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Preface 


"It  is  not  possible  to  predict  with  complete 
accuracy  and  confidence  the  expected  impact 
area  of  an  artificial  Earth  satellite  after 
unaided  orbit  decay.  In  fact,  it  has 
usually  been  impossible  to  ascertain  the 
locations  of  the  remains  of  most  satellites 
after  their  descent  from  orbit.”  (1) 

This  quotation  is  essentially  as  valid  today  as  it  was  in  1965  when  it 
appeared  in  the  report  of  a  joint  study  by  the  USAF  SPACETRACK  System, 
the  Smithsonian  Astrophysical  Observatory,  and  the  RAND  Corporation. 

This  is  particularly  true  when  one  considers  the  actual  limitations 
associated  with  the  term  "unaided".  An  unaided  orbital  decay  repre¬ 
sents  a  natural  decay  from  Earth  orbit  due  to  the  action  of  atmospheric 
drag  in  a  very  ill-defined  dynamics  environment.  The  uncertainties 
associated  with  this  process  include  large,  time  and  space  dependent, 
variations  in  atmospheric  density  and  violent  changes  in  vehicle 
configuration  and  structure  due  to  the  ablative  and  dynamically 
unstable  process  of  reentry.  With  many  hundreds  or  thousands  of 
artificial  Earth  satellites  and  pieces  of  launch  booster  debris 
currently  in  orbit  at  various  orbital  inclinations,  it  is  extremely 
difficult  to  determine  and  predict  the  reentry  impact  locations. 

The  underlying  limitation  of  determining  the  reentry  trajectory 
resides  within  the  intractable  mathematical  nature  of  the  true 
dynamic  processes  of  reentry.  Standard  atmosphere  models  often  provide 
valuable  information  on  the  magnitude  of  the  atmospheric  density  and 
its  functional  variation  with  altitude.  Unfortunately,  these  models  are 
developed  to  represent  a  "mean"  model  of  the  atmosphere  and,  in  the 
reentry  altitude  regions  are  usually  limited  to  data  collected  in  the 
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northern  hemisphere.  As  such,  the  mean  models  of  a  standard  atmosphere 
do  not  properly  represent  the  density  variations  along  each  individual 
satellite  decay  trajectory.  The  decaying  reentry  body  is  also  subject 
to  a  wide  variation  in  geometric  forms  as  it  passes  through  the 
atmosphere.  Factors  such  as  ablation,  random  rotation,  and  fragmenta¬ 
tion  preclude  definition  of  the  vehicle  configuration  throughout  the 
reentry . 

These  circumstances  present  difficult  problems  if  one  attempts 
to  determine  the  reentry  trajectory  by  developing  a  dynamic  force 
expression  from  standard  aerodynamic  theory.  An  accurate  solution  to 
the  hypersonic  flow  field  surrounding  the  reentry  body  is  not  possible 
in  an  environment  where  one  cannot  specify  the  atmospheric  properties, 
the  time  varying  character  of  the  vehicle  geometry,  and  the  location 
and  shape  of  the  reentry  shock  wave. 

The  complexities  associated  with  these  processes  force  one  to 
consider  using  a  more  simplified  representation  of  the  reentry  dynamics. 
The  inadequacies  of  these  simple  modeling  techniques  are  evident  in  the 
current  operational  techniques  of  propagating  the  final  orbit 
determination  state  vector  and  state  covariance  matrix  through  to  Earth 
impact.  Uncertainties  of  many  hundreds  or  thousands  of  nautical  miles 
along  the  trajectory  ground  trace  are  apparent  in  the  ultimate  impact 
location. (1) 

Assuming  the  availability  of  trajectory  observations  during 
reentry  from  an  orbital  sensor,  one  essentially  has  a  global  visibility 
of  the  arbitrary  decay  trajectories.  With  a  simplified  definition  of 
the  physical  dynamics  of  reentry,  one  can  incorporate  the  observations 
into  an  estimation  scheme  to  improve  the  knowledge  of  the  ultimate 
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Earth  impact  point.  With  frequent  enough  observations,  one  can 


potentially  use  a  model  which  would  be  inadequate  over  a  long  period  of 
time.  The  principal  goal  of  this  research  was  the  development  and 
application  of  such  an  estimation  technique.  A  linearized,  differ¬ 
ential  corrector  method  is  utilized  which  is  directly  applicable  to 
operational  agencies  such  as  the  USAF  SPACETRACK  System  which  must 
reconstruct  the  impact  locations  of  many  reentry  objects  with  minimum 
a  priori  knowledge  of  the  given  reentry  dynamics. 

The  research  approach  was  predicated  on  extending  the  existing 
orbit  determination  methods  into  the  reentry  regions.  This  included 
the  use  of  specified  observations  from  an  orbital  sensor  assumed  to 
provide  angular  (azimuth  and  elevation)  measurements  from  an  infrared 
(IR)  source  at  a  fixed,  discrete  data  rate.  The  major  limitations  in 
direct  application  of  a  differential  corrector  with  a  deterministic 
dynamics  model  will  be  shown  to  be  variations  and  uncertainties 
associated  with  the  reentry  dynamics.  Standard  estimation  methods 
would  consider  adding  a  pseudo-noise  compensation  to  the  dynamic 
equations  of  motion  to  incorporate  the  dynamic  uncertainties  present. 
Unfortunately,  the  operational  experience  of  agencies  such  as  the  USAF 
SPACETRACK  System  shows  that  up  to  100  orbital  revolutions  of 
empirical  tracking  data  have  been  required  for  proper  "tuning"  of  such 
a  pseudo-noise  compensation  technique.  The  modified  dynamics  have  then 
been  found  to  apply  only  to  the  given  satellite  in  question. 

In  the  more  uncertain  dynamics  regions  of  reentry,  with  a 
single  short  arc  of  empirical  data  from  trajectory  observations,  even 
less  potential  exists  to  apply  a  pseudo-noise  compensation  successfully 
to  the  dynamics  model  of  a  given  reentry.  A  basic  contribution  of  this 


research  will  be  shown  to  be  the  application  of  an  "ad  hoc"  scalar 
fading  memory  parameter  which  is  adaptively  determined  from  the 
estimation  process.  This  represents  an  extension  of  the  earlier  work 
of  Morrison  (2)  and  Sorrenson  and  Sacks  (3)  to  the  reentry  estimation 
application.  This  scalar  parameter  is  used  to  multiply  tne  terms  of 
the  estimator-computed  state  covariance  matrix  prior  to  an  observation 
update.  This  new  matrix  is  referred  to  as  a  "deweighted"  state 
covariance  matrix.  The  scalar  deweighting  parameter  provides  a  fading 
memory  on  the  previous  history  of  observations. 

The  magnitude  of  the  fading  memory  parameter  is  selected  to 
force  acceptable  estimator  performance  relative  to  this  "deweighted" 
state  covariance  matrix.  This  acceptable  performance  is  demonstrated 
via  a  Monte  Carlo  analysis.  The  mean  error  in  the  estimator  solutions 
(mean  state  estimate  -  true  state  solution)  is  much  smaller  than  the 
standard  deviations  from  the  estimator-computed  deweighted  state 
covariance  matrix.  The  RSS  (root  sum  square)  of  the  mean  square  posi¬ 
tion  and  velocity  errors  in  the  Monte  Carlo  solutions  compares  closely 
to  the  standard  deviations  from  the  deweighted  state  covariance  matrix. 

With  its  successful  application  to  reentry  trajectories,  the 
proposed  technique  provides  a  valuable  tool  for  astrodynamic  research. 
In  the  near  term,  estimated  Earth  impacts  with  reasonably  valid 
uncertainties  can  quickly  be  computed  after  orbit  decay.  In  the  long 
term,  the  method  offers  an  ability  to  construct  an  empirical  data  base 
of  reentry  trajectories  from  which  a  pseudo-noise,  or  adaptively 
determined  pseudo-noise,  compensation  method  can  potentially  be  devel¬ 
oped.  This  should  afford  a  better  means  to  understand  the  true 


dynamics  of  the  decay  trajectories,  such  that  more  sophisticated  estima 
tion  techniques  can  successfully  be  applied  to  this  problem. 
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Abstract 

A  technique  was  developed  for  estimation  of  decayed  Earth  satel¬ 
lite  reentry  trajectories  to  provide  impact  locations  for  debris  search 
requirements.  The  technique  used  a  linearized  differential-corrector 
as  an  extension  of  existing  orbit  estimation  methods.  The  reentry 
observations  consisted  of  angular,  infrared  measurements  from  orbital 
sensors.  An  eight  dimension  state  vector  was  used  with  components  for 
position,  velocity,  a  ballistic  parameter,  and  the  scale  height  from 
an  isothermal  density  model.  Simulated  data  runs  identified  the 
uncertain  dynamics  of  the  true  reentry  process  as  the  most  significant 
impact  on  estimator  performance.  The  uncertain  dynamics  pose  signifi¬ 
cant  problems  for  standard  model  compensation  methods  such  as  adaptive 
pseudo-noise  compensation  or  more  sophisticated  techniques  such  as 
statistical  linearization  or  higher  order  filters.  The  adaptive 
determination  of  an  "ad  hoc"  scalar  fading  memory  parameter  was  used  to 
modify  the  estimator-computed  state  covariance  matrix.  The  bias  in 
the  state  estimates  were  well  within  the  variance  from  this  modified 
covariance  matrix.  The  standard  deviations  from  the  modified 
covariance  compare  closely  to  the  root-mean-square  errors  of  the  true 
solution  over  a  range  of  simulated  truth  model  data.  The  uncertainties 
associated  with  the  impact  locations  are  anticipated  to  be  at  least  one 
order  of  magnitude  better  than  propagation  of  the  final  orbit  estima¬ 
tion  covariance  matrix  to  Earth  impact. 
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Chapter  I  -  Introduction 


A.  Research  Overview 

The  objective  of  this  research  was  to  develop  an  application  of 
an  estimation  technique  which  would  support  debris  search  requirements 
in  the  impact  area  of  decayed  Earth  satellites.  This  necessitated 
determining  impact  locations  and  defining  the  uncertainties  associated 
with  those  positions.  The  engineering  considerations  included  the  use 
of  specified  orbital  observer(s)  providing  angular  infrared  (IR)  tra¬ 
jectory  observations  at  fixed,  discrete  time  intervals  and  large 
uncertainties  in  the  mathematical  character  of  the  true  reentry 
dynamics.  Angular  observations  were  assumed  since  they  could  be 
obtained  from  orbital  satellites  which  view  the  reentry  using  a  passive 
IR  sensor.  The  aerodynamic  heating  of  the  reentry  satellite  provides 
a  convenient  infrared  signature  as  it  decays  through  the  Earth  atmos¬ 
phere.  While  the  additional  use  of  range,  or  range-rate,  observations 
may  improve  the  observability  of  the  reentry  trajectory,  these  would 
necessitate  a  more  complex  active  sensor  system,  such  as  a  radar.  The 
power  required  for  radar  tracking  varies  as  a  function  of  range  to  the 
fourth  power  between  the  radar  and  the  reentry  satellite.  These 
sensors  would  be  much  more  costly  to  deploy  In  orbit  to  provide  full 
visibility  of  the  Earth. 

The  uncertain  dynamics  of  the  reentry  process  present  many 
difficulties  for  existing  estimation  techniques.  A  significant  portion 
of  this  research  was  to  identify  the  limits  of  the  various  estimation 
techniques  which  might  be  considered  for  application  to  this  problem. 
The  se' acted  approach,  which  uses  an  adaptively  determined,  ad  hoc. 
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scalar  fading  memory  parameter,  is  offered  as  an  interim  solution  to 
the  problem.  The  research  approach  concentrated  on  extending  existing 
orbit  determination  methods  into  the  reentry  region.  These  techniques 
generally  use  a  deterministic  dynamics  model  and  linearize  about  a 
reference  trajectory  in  a  Taylor's  series  expansion.  The  use  of  a 
linearized  approach  was  retained  in  this  application  to  maintain 
consistency  with  the  existing  orbit  estimation  methods.  The  line¬ 
arized  approach  also  allows  one  to  concentrate  on  improving  the  first 
order  effects  of  the  uncertain  reentry  dynamics.  Once  these  are  better 
understood,  a  more  accurate  dynamics  model  may  potentially  be  defined. 
This  more  accurate  model  may  require  application  of  a  more  sophisti¬ 
cated  estimation  technique  to  handle  the  higher  order  effects  properly 
within  the  reentry  dynamics.  Methods  such  as  statistical  linearization 
and  nonlinear  estimation  are  best  suited  to  applications  where  a  good 
definition  exists  for  the  functional  form  of  the  nonlinear  terms  in 
the  dynamics  model. 
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B.  Background 


Various  operational  space  agencies  -  the  National  Aeronautics 
and  Space  Administration  (NASA),  the  USAF  Space  Defense  Center 
(SPACETRACK) ,  and  the  USAF  Satellite  Control  Facility  (AFSCF)  -  have  a 
mission  responsibility  to  determine  impact  locations  of  decayed  Earth 
satellites-  Most  of  their  respective  orbit  computation  systems,  with 
minor  differences,  have  a  similar  kind  of  data  processing  capability. 

A  weighted  least-squares,  differential  correction  algorithm  is  used  to 
process  observational  data  to  calculate  an  "optimal"  estimate  of  the 
satellite  orbit.  In  the  appropriate  altitude  regimes,  the  estimation 
algorithm  is  capable  of  computing  an  estimate,  for  both  the  state 
trajectory  positions  and  velocities,  and  various  perturbations  to  the 
basic  orbital  motion  (i.e.,  ballistic  parameters,  perturbations  to 
mean  atmospheric  density  values,  etc.)  (4,5,6).  Applications  of  the 
differential  correction  technique  have  evolved  from  astrodynamic 
applications  where  the  basic  assumptions  of  the  methods  apply.  These 
include:  i)  a  deterministic  dynamics  model,  ii)  noise  corrupted 
observation  measurements,  and  iii)  batch  processing  of  a  large  number 
of  obser  ations  at  selected  points  in  time  to  develop  corrections  and 
covariance  values  for  the  state  variables  at  an  arbitrarily  selected 
point  along  the  reentry  trajectory  referred  to  as  an  "epoch".  The 
estimates  of  the  state  vector  and  covariance  are  propagated  forward  to 
a  new  epoch  point  where  new  observations  are  batch  processed  to  update 
the  trajectory  estimation. 

In  the  orbital  regime,  this  technique  has  met  with  considerable 
success  due  to  a  number  of  factors: 
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1)  The  availability  of  sufficient  and  adequately  spaced 
observations  from  ground  based  radar  or  telemetry 
measurements . 

2)  Relatively  well  defined  satellite  geometric  configurations 
and  relative  insensitivity  to  the  small  aerodynamic  effects 
of  a  rarefied  atmosphere  which  allows  one  to  use  a 
deterministic  dynamics  model. 

3)  The  basic  validity  of  linearizing  about  a  reference 
trajectory  while  processing  observations  in  a  batch 
processing  mode  (where  all  observations  are  processed 
together  for  an  orbit  update)  or  a  sequential  batch 
processing  mode  (which  processes  a  subset  of  the 
observations  of  various  sizes,  along  the  trajectory  or 
orbit) . 

As  one  enters  the  reentry  region  (after  a  spiral  decay  from 
orbit  has  commenced  and  at  altitudes  generally  well  below  125  KM),  the 
ability  to  extend  the  existing  orbit  estimation  techniques  becomes 
appreciably  more  difficult.  This  difficulty  is  a  direct  consequence 
of  the  following  factors: 

1)  The  lack  of  universally  available  tracking  observations 
of  the  final  decay  trajectory, 

2)  Significant  uncertainties  in  the  vehicle  dynamics  and  in 
the  atmospheric  density,  which  limit  the  use  of  a 
deterministic  dynamics  model,  and 

3)  A  significantly  more  nonlinear  set  of  dynamics  within  the 
equations  of  motion. 


Consider  the  aerodynamic  acceleration  expression  most  frequently 
used  in  the  orbital  equations  of  motion.  The  dominant  aerodynamic 
acceleration  term  is  due  to  atmospheric  "drag”  acting  along  the 
instantaneous  velocity  vector  of  the  satellite  motion: 


^  2  P  VR  S 


where : 


a^  =  Acceleration  due  to  "drag",  generally  defined  to  be 

the  resolution  of  the  aerodynamic  accelerations  along 
the  three  spatial  coordinates  into  one  vector,  opposite 
in  direction  to,  and  along  the  instantaneous  velocity 


p  =  Atmospheric  density 

V  =  Scalar  magnitude  of  the  velocity  vector  relative  to  the 
K 

rotating  atmosphere 

V  =  A  column  vector  of  the  three  components  of  the  velocity 
K 

relative  to  the  rotating  atmosphere 
6  =  "Ballistic  parameter"  which  equals  C^A/M,  with 

Cp  =  "Drag"  coefficient 
A  =  Satellite  reference  area 
M  =  Satellite  mass 

In  the  reentry  regions,  the  uncertainties  associated  with  the 
geopotential  forces  are  orders  of  magnitude  smaller  than  the 
uncertainties  of  the  aerodynamic  forces.  In  a  relative  sense,  the 
non-aerodynamic  forces  can  be  considered  as  deterministic  geopotential 
forces.  These  are  quite  adequately  derived  from  a  well  developed  set 
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of  zonal  and  tesseral  harmonics  associated  with  the  geopotential 
function  of  a  standard  Earth  model  (7). 


The  two  basic  quantities  which  induce  uncertainty  into  the  aero¬ 
dynamic  acceleration  expression  are  8  and  p.  Because  these  two  quanti¬ 
ties  appear  as  the  product  8p  within  the  dominant  acceleration 
expression,  the  result  of  uncertainty  in  either  is  magnified  in  this 
multiplicative  fashion. 

As  Reference  1  reports,  a  change  of  only  ±10%  in  either  8  or  p 
at  the  point  of  orbital  decircularization,  and  the  onset  of  spiral 
decay  from  orbit,  can  ultimately  result  in  up  to  a  one-half  Earth 
revolution  difference  in  the  ultimate  impact  location  of  the  reentry 
satellite,  when  not  updated  by  further  observations  beyond  this  reentry 
point.  This  fact  is  particularly  apparent  if  one  thinks  of  the 
sensitivity  of  the  direction  and  magnitude  of  the  velocity  vector  at 
the  onset  of  decay  to  variations  in  either  8  or  p .  These  uncertain¬ 
ties  are  direct  contributors  to  the  many  hundreds  or  thousands  of 
nautical  mile  uncertainties  of  the  impact  location  when  obtained  by  a 
straightforward  propagation  of  the  last  orbital  vector  and  covariance. 
B.l.  Trajectory  Observations: 

With  the  availability  of  sufficiently  dense  and  accurate 
observations,  significant  improvements  to  the  decayed  satellite 
trajectory  estimation  should  be  available.  Tracking  data  from  Earth 
surface  sites  is  generally  limited  due  to:  i)  wide  geometric 
separation  of  a  very  few  Earth  tracking  stations,  and  ii)  short 
trajectory  spans  of  observability  yielding  few  observations  at  any 
site(s)  which  may  fortunately  be  positioned  to  view  the  reentry.  If 
one  postulates  a  network  of  synchronous  tracking  satellites,  orbitally 


positioned  to  provide  a  global  visibility  of  the  Earth,  the  limited 
availability  of  observations  for  the  general  decay  trajectory  from 
ground  sources  should  be  eliminated  or  minimized. 

This  orbitally  positioned  sensor  is  assumed  to  provide  angle 
only  IR  observations  (i.e.,  no  range  measurement)  at  a  fixed  data  rate 
of  every  10  seconds.  The  implications  of  these  assumptions  will  be 
examined  in  Chapters  II  and  III.  With  the  availability  of  such  obser¬ 
vations  in  azimuth  and  elevation,  one  may  define  a  coordinate 
reference  frame  as  shown  in  Figure  1. 

B.2.  Dynamic  Uncertainties: 

B.2.1.  Vehicle  Uncertainties: 

The  uncertainties  in  the  ballistic  parameter,  6,  result  from  a 
variety  of  sources  for  the  arbitrary  orbital  decay.  These  include 
unknown  spacecraft  orientation,  asymmetric  vehicle  ablation,  and 
vehicle  fragmentation.  In  the  orbital  regime,  ablation  and  fragmenta¬ 
tion  are  not  factors.  In  the  reentry  regime,  high  deceleration  and 
aerodynamic  heating  occur.  Abrupt  and  rapidly  changing  variations 
result  due  to  ablation  and  structural  fragmentation.  These  effects 
are  further  complicated  by  their  coupling  with  the  atmospheric  density 
term,  p,  and  their  irregular  occurrence  along  the  reentry  trajectory. 

In  the  orbital  regime  there  usually  exists  a  well  defined  space¬ 
craft  angular  momentum  vector  derived  from  the  spacecraft  design  or 
developed  from  observations  of  the  satellite  or  booster  debris.  This 
allows  definition  of  at  least  a  "mean"  value  of  the  ballistic 
parameter,  6,  for  specified  segments  of  the  orbit.  Even  without 
fragmentation,  the  behavior  and  definition  of  the  angular  momentum 
vector  is  difficult  to  model  in  the  highly  variable  force  field 
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experienced  during  reentry. 

B.2.2.  Atmospheric  Density  Uncertainties: 

The  uncertainties  in  atmospheric  density  are  due  to  several 
well  documented  perturbations  to  the  mean  values  of  the  standard  den¬ 
sity  models  (8).  These  include: 

1)  Seasonal  variations 

2)  Diurnal  variations 

3)  Solar  radiation  effects 

4)  Magnetic  storm  effects 

5)  Auroral  effects 

6)  Local  climatology  effects 

The  magnitude  of  these  departures  from  the  mean  density  pro¬ 
files  is  well  documented  in  Reference  9,  which  shows  ±95%  certainty 
deviations  of  up  to  ±60%  from  the  mean  density  of  the  standard 
atmosphere. 

The  present  understanding  of  the  mathematical  structure  and 
underlying  physics  of  the  standard  atmosphere  models  is  well  docu¬ 
mented  in  a  number  of  references  (9,10,11,12).  However,  the  effects 
of  these  variations  are  both  difficult  to  quantify  and  difficult  to 
model  on  a  spatial  and  temporal  basis  for  the  general  decay  trajectory. 
This  is  a  consequence  of  a  number  of  factors.  There  generally  is  a 
lack  of  detailed  density  measurements  along  the  individual  reentry 
trajectory.  These  could  be  used  to  establish  boundary  conditions  for 
the  mathematics  which  describe  the  propagation  of  the  perturbations 
through  the  atmosphere.  Most  mean  density  models  have  been  developed 
from  a  limited  number  of  geographically  distributed  density  measure¬ 
ments.  Most  measurements  below  orbital  altitudes  and  above  normal 


aircraft  flight  altitudes  are  from  sounding  rocket  or  balloon 
instrumented  flights  densely  grouped  near  45°  North  latitude.  The 
mean  density  models  are  more  accurate  in  this  region.  Accurate 
specification  of  density  is  required  for  the  particular  reentry  trajec¬ 
tory  being  estimated,  and  this  reentry  may  occur  anywhere  around  the 
Earth . 

While  many  operational  space  agencies  utilize  corrections  to  the 
standard  atmosphere  density  values,  these  corrections  are  difficult  to 
determine  due  to  their  wide  variation  over  time  in  the  reentry 
altitude  regions.  In  the  orbital  applications,  the  three  density 
variations  generally  added  to  correct  the  mean  density  values  include 
(4,5,6): 

1)  Those  due  to  magnetic  flux  variations,  Ao 

m 

2)  Those  due  to  solar  radiation  variations,  Ac 

s 

3)  Those  related  to  the  geometry  relationship  between  the 
satellite  position  and  the  "diurnal  bulge"  of  the  Earth 
atmosphere . 

This  area  of  most  severe  solar  heating  occurs  along  the  vector  between 
the  sun  and  Earth  center.  The  intersection  of  this;  vector  with  the 
Earth  surface  is  referred  to  as  the  solar  subpoint,  or  nadir. 

In  the  orbital  applications,  the  first  two  modifications  are 
usually  added  to  the  standard  density  values  as  a  constant  over  a 
segment  of  the  orbit.  The  magnitude  of  these  modifications  are 
determined  from  one  of  two  separate  computations.  If  they  exist  in  the 
orbital  altitudes  and  inclinations  of  interest,  special  "calibration" 
satellites  (whose  aerodynamic  configurations  and  orbits  are  precisely 
known)  are  used  to  determine  the  magnitude  of  the  density  modifications. 
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Alternately  (often  using  the  values  from  the  calibration  technique  as 
initial  conditions),  the  density  modifications  are  added  to  the  state 
vector  and  separately  estimated  along  with  the  orbital  trajectory  of 
interest.  In  either  case  this  increases  the  complexity  of  the  equations 
of  motion.  If  the  modifications  are  also  estimated  along  with  the 
satellite  trajectory,  this  increases  the  dimension  of  the  state  vector. 
For  global  applications  in  the  reentry  region,  there  is  no  systematic 
means  to  establish  the  initial  conditions  of  these  density  modifica¬ 
tions  which  are  valid  over  all  altitude  regions  of  the  general  decay 
traj  ectory . 

By  default,  therefore,  attempts  to  apply  the  orbital  methods  to 
the  reentry  regime  have  often  been  confined  to  estimating  (3  along  with 
the  position  and  velocity  variables.  A  state  vector  of  seven 
dimensions  is  used,  with  atmospheric  density  derived  from  a  standard 
density  model.  This  approach  attempts  to  group  all  the  unmodeled 
parameters  within  3.  This  approach  is  limited  by  attempting  to  estimate 
a  variety  of  unknown  functional  forms  as  a  constant  over  some  short 
trajectory  span.  This  motivates  one  to  consider  a  finite  memory  or 
fading  memory  estimator.  However,  attempts  to  estimate  many  additional 
physical  parameters  may  also  yield  identification  or  observability 
problems,  particularly  when  observations  are  limited  to  angular 
measurements  only.  The  implications  of  these  limitations  will  be 
addressed  in  Chapter  II  with  the  proposed  estimator  formulation. 

With  the  above  considerations  in  mind,  an  eight  dimensional  state 
vector  was  chosen  for  the  reentry  application.  The  basic  three 
position  and  three  velocity  components  were  augmented  by  the  ballistic 
parameter,  8,  and  the  density  scale  height,  Q,  from  an  exponential 


density  model.  The  azimuth  and  elevation  measurements  from  the 
orbiting  sensor(s)  were  then  incorporated  within  a  differential 
corrector  to  estimate  the  reentry  trajectory.  The  discussion  below 
will  summarize  the  research  which  is  presented  in  more  detail  in  the 
subsequent  chapters. 


C.  Research  Approach 


This  report  will  first  document  the  model  definition  consid¬ 
erations  and  apply  the  basic  differential  corrector  formulation  to  a 
variety  of  simulated  observations.  Modifications  to  the  existing 
orbit  determination  methods  are  related  to  the  basic  underlying 
principles  of  estimation  theory  throughout  the  text.  One  must 
consider  the  limitations  of  operating  in  a  regime  of  uncertain  dynamics 
with  an  inability  to  exploit  the  flexibility  of  many  standard  model 
compensation  methods  of  estimation  theory.  The  physical  circumstances 
and  factors  which  present  these  limitations  are  discussed,  as 
appropriate,  within  the  subsequent  chapters. 

Using  the  ballistic  parameter  and  the  density  scale  height  of 
an  exponential  density  model  to  augment  the  position  and  velocity 
variables,  a  linearized  technique  was  applied  to  estimate  the  reentry 
trajectories.  The  derivation  of  the  basic  linearized  differential 
corrector  and  reentry  dynamics  model  is  provided  in  Chapter  II. 

Careful  definition  of  the  estimator  dynamics  model  was  neces¬ 
sary  to  insure  a  valid  linearization  relative  to  the  reference  reentry 
trajectory  and  valid  specification  of  initial  conditions  for  the  state 
variables  and  state  covariance  matrix  entries.  Even  with  a  simplified 
dynamics  model,  the  differential  corrector  technique  must  be  examined 
to  insure  valid  linearization  for  both  state  and  covariance  updates 
and  state  covariance  propagation  over  time.  Restrictions  on  the  time 
span  of  valid  linearization  over  the  trajectory  have  also  been 
examined  as  part  of  this  research.  The  ultimate  effect  of  limiting 
this  time  span  is  to  restrict  the  number  of  observations  one  may 
incorporate  into  a  single  update  at  the  reference  time  epoch.  In 


regions  of  high  dynamic  uncertainty,  the  linearization  assumptions  of 
the  deterministic  trajectory  dynamics  may  require  use  of  a  sequential 
processor  (single  observation  updates).  The  implementation  of  a 
linearization  validity  check  is  presented  in  Chapter  II. 

The  physical  environment  within  which  the  estimator  would 
ultimately  operate  was  also  given  special  attention.  Extensive 
simulation  analyses  were  completed  to  examine  estimator  performance 
when  subjected  to  simulated  observations  derived  from  a  realistic 
"truth  model".  These  analyses  examined  the  effects  of  mismatch 
between  the  truth  model  and  the  estimator  model  dynamics,  accuracy 
variations  on  the  angular  observations,  multiple  orbital  observation 
locations,  and  variations  of  the  geometric  rela^ic.-thips  between  the 
observing  satellite(s)  and  the  reentry  trajectory.  The  truth  model 
selections  were  chosen  to  simulate  a  full  range  of  circumstances  which 
the  estimator  application  would  ultimately  encounter. 

A  series  of  Monte  Carlo  analyses  are  presented  which  identify 
the  model  dynamics  as  the  most  important  item  impacting  the  estimator 
performance.  This  assessment  was  accomplished  by  first  considering  the 
mean  bias  of  the  trajectory  position  and  velocity  estimates  relative  to 
the  magnitude  of  the  standard  deviations  from  the  estimator-computed 
state  covariance  matrix.  A  comparison  was  also  made  between  the 
magnitudes  of  the  standard  deviations  from  the  estimator-computed  state 
covariance  matrix  and  those  derived  from  the  Monte  Carlo  samples.  This 
allows  an  examination  of  the  systematic  error  in  the  state  estimate  and 
the  validity  of  the  random  error  indicated  by  the  estimator  state 


covariance  matrix. 


With  matched  dynamics  in  the  estimator  and  the  truth  model,  the 


t$ias  in  the  estimator  solution  remains  small  relative  to  the  standard 
deviations  from  the  estimator-computed  state  covariance  matrix.  Also, 

i 

j  the  standard  deviations  from  the  estimator-computed  covariance  compare 

closely  to  those  derived  from  the  Monte  Carlo  data.  With  a  mismatch 

j  between  the  estimator  and  the  truth  model  dynamics,  the  bias  and  the 

ij 

!  standard  deviations  derived  from  the  Monte  Carlo  analyses  grow  large 

relative  to  the  standard  deviations  from  the  estimator-computed  state 
covariance  matrix. 

| 

I  A  review  of  several  model  compensation  methods  is  presented 

relative  to  this  application.  These  include  techniques  for  adding  a 
pseudo-noise  compensation  to  the  model  dynamics,  adaptive  estimation 

i 

methods,  and  state  covariance  deweighting  techniques.  The  limitations 
associated  with  the  use  of  each  of  these  methods  for  the  current 

f 

|  reentry  application  are  discussed. 

I 

Chapter  III  presents  the  derivation  of  a  fading  memory 
differential  corrector  analogous  to  the  linear  estimation  developments 
of  previous  researchers  (2,3).  Applications  are  made  to  a  wide  class 
of  dynamic  variations  in  the  reentry  trajectory.  Successful  applica¬ 
tion  is  achieved  by  an  adaptive  determination  of  an  ad  hoc  scalar  mul¬ 
tiplier,  y.  This  scalar  parameter  is  used  to  multiply  the  terms  of  the 
state  covariance  matrix  prior  to  an  observation  update,  and  thus, 
implements  a  fading  memory  on  the  processing  of  earlier  observations. 
This  is  accomplished  by  examination  of  the  size  of  the  change  in  state 
variables,  <5x^,  at  each  update  epoch  along  the  reentry  trajectory. 

The  magnitudes  of  each  6x^  are  compared  to  the  magnitude  of  the 
standard  deviation  of  their  respective  terms  in  the  deweighted  state 
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covariance  matrix,  as  computed  in  the  estimation  algorithm  itself. 

This  affords  a  measure  to  determine  which  trajectory  time  points 
require  application  of  fading  memory  within  the  estimator. 

A  series  of  Monte  Carlo  results  are  presented,  showing  the 
application  of  this  adaptively  determined,  «.  I  hoc  scalar  fading  memory 
parameter  for  estimation  of  an  anticipated  true  reentry  dynamics 
trajectory.  An  examination  of  bias  within  tht  estimator  solution  with 
these  simulated  true  dynamics  is  provided.  The  limitations  associated 
with  using  the  estimator-computed  covariance  matrix  for  impact 
uncertainties  are  also  discussed. 

Finally,  a  tangent  plane  projection  of  the  Earth  impacts  of 
several  estimator  applications  will  be  shown,  demonstrating  the 
ability  of  the  fading  memory  method  to  provide  impact  locations  with 
bias  magnitudes  within  the  standard  deviations  of  the  "deweighted" 
covariance  matrix.  The  standard  deviation  of  the  position  error  from 
the  deweighted  state  covariance  matrix  provides  a  good  definition  of 
the  uncertainties  in  the  estimated  Earth  impact  location,  thus  it  can 
be  used  to  define  a  search  area  for  the  recovery  of  the  satellite 
debris.  These  results  illustrate  the  viability  of  the  method  to 
estimate  decayed  satellite  impact  locations  and  uncertainties 
significantly  improved  over  existing  astrodynamic  applications. 

The  basic  contribution  of  the  research  is  to  provide  a 
technique  which  will  help  to  determine  the  satellite  impact  locations 
with  reasonable  specification  of  their  uncertainties  in  an  Earth 
tangent  plane  coordinate  system  projection.  Standard  estimation 
techniques  are  applied,  but  in  a  unique  manner  as  dictated  by  the  needs 
of  this  specific,  difficult  and  heretofore  inadequately  resolved 
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problem  area.  Important  highlights  of  the  research  include  the 
identification  of  the  limitations  associated  with  using  the  standard 
differential  corrector  for  this  application.  The  research  also 
documents  the  problems  which  the  decayed  reentry  problem  poses  for 
more  sophisticated  estimation  techniques  such  as  pseudo-noise 
compensation  or  statistical  linearization. 

A  secondary  goal  is  achieved  by  providing  an  estimation 
approach  which  can  be  used  to  develop  a  large  empirical  data  base  of 
decayed  trajectories.  These  may  potentially  be  used  to  develop  the 
alternative  of  a  pseudo-noise  dynamic  compensation  to  the  model 
dynamics.  The  estimator  developed  represents  the  first  successful 
exploitation  of  global  tracking  observations  from  an  orbital  sensor, 
while  simultaneously  estimating  an  update  of  the  atmospheric  density 
with  angular  observations. 


Chapter  11  -  The  Reentry  Estimator 


As  mentioned  in  the  Introduction,  the  initial  research  approach 
concentrated  on  extending  the  existing  orbit  estimation  techniques  to 
the  decayed  satellite  problem.  These  methods  generally  use  a 
deterministic  dynamics  model  within  a  differential  corrector.  This 
chapter  will  document  the  successes  and  limitations  of  using  such  a 
technique,  considering  the  dynamics  uncertainties  present  in  this 
application. 

The  first  step  in  applying  the  differential  corrector  to  the 
general  decay  trajectory  estimation  problem  involves  definition  of 
the  reentry  dynamics  model  and  development  of  the  estimator 
structure.  The  dynamics  model  considerations  were  introduced  in 
Chapter  I.  This  chapter  will  begin  with  a  derivation  of  the  basic 
differential  corrector  as  it  will  be  applied  to  the  estimation  of  the 
reentry  trajectory  and  Earth  impact.  A  complete  definition  of  the 
dynamics  model  is  also  presented.  This  includes  geopotential 
accelerations  from  the  Smithsonian  Astrophysical  Observatory  SAO-III 
Earth  Model  (7) ,  utilization  of  an  exponential  atmosphere  density 
based  upon  a  least-squares  fit  to  the  U.S.  1962  Standard  Atmosphere, 
and  the  assumption  of  a  constant  ballistic  parameter,  8,  over  a  given 
segment  of  the  reentry  trajectory.  The  implementation  details  of  the 
estimator  will  then  be  discussed. 

Numerical  simulation  results  which  help  establish  the  limits 
of  the  estimator  performance  are  presented.  These  include  single 
sample  simulation  runs,  selected  Monte  Carlo  simulation  results,  and 
special  r-iraerical  examples  peculiar  to  the  observer  geometry 
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relationships  and  reentry  dynamics  of  this  application.  The  overall 
impact  of  the  Chapter  II  results  will  be  to  show  the  basic  differential 
corrector  with  deterministic  dynamics  as  inadequate  for  processing  true 
reentry  data.  The  deterministic  dynamics  and  infinite  memory  of  the 
basic  formulation  cause  the  estimator  to  yield  significantly  biased 
solutions  relative  to  the  standard  deviations  of  the  estimator- 
computed  state  covariance  matrix.  These  occur  when  processing  data 
reflecting  the  dynamic  variations  anticipated  in  true  reentry 
trajectories.  The  use  of  the  estimator  in  this  form  would  not  provide 
an  accurate  estimate  of  Earth  impacts  for  satellite  debris  search 

requirements.  With  a  large  systematic  error  (bias)  in  the  state  i 

1 

estimate,  the  estimator  state  error  covariance  matrix  could  not  be  j 

| 

utilized  to  define  a  search  area  around  the  impact  location.  The  j 

position  bias  from  the  estimator  solution  may  be  significantly  removed  j 

from  the  actual  impact  location.  It  may  also  be  much  larger  in 

magnitude  than  the  random  error  indicated  by  the  estimator  state 

covariance  matrix.  The  results  of  Chapter  III  will  address  this 

limitation. 


A.  The  Weighted  Least-Squares  Differential  Corrector 

The  basic  differential  corrector  includes  a  nonlinear  dynamics 

model : 

x  =  f (x,t)  (2) 

where  f(x,t)  is  a  deterministic  function  of  the  state  variables  and  a 
continuous  function  of  time.  The  overbar,  x,  denotes  a  vector 
quantity.  Such  a  deterministic  dynamics  model  is  appropriate  when 
f(x,t)  provides  a  good  definition  of  the  true  dynamics. 

A  set  of  discrete  observations  are  related  to  the  state 
variables  by  the  general  nonlinear  relationship: 

z(t  )  =  h(x(t  ),t  )  +  r (t  )  (3) 

n  n  n  n 

where  r(t  )  is  a  zero  mean,  corruptive  noise  term  of  covariance  R(t  ) 
n  n 

at  the  observation  time  t  . 

n 

Assume  the  availability  of  a  reference  trajectory,  x^Ct),  with 

initial  conditions  x  (t  ) ,  at  the  epoch  time,  t  .  One  seeks  to  find  a 

o  o  o 

correction,  6x(t)  ,  which  will  minimize  a  weighted  quadratic  cost  func¬ 
tion  of  the  observation  residuals. 

Assume  that  the  true  dynamics  solution  can  be  represented  by: 

x(t)  =  XQ(t)  +  6x(t)  (4) 

The  reference  trajectory,  ,  satisfies  the  dynamics  model  as: 

x  (t)  =  f  (x  (t)  ,t)  (5) 

o  o 

under  the  premise  that  the  6x(t)  is  a  "small"  deviation  from  the  true 


solution,  x(t).  If  one  expands  Equation  2  in  a  Taylor’s  series  about 
the  reference  trajectory,  the  following  results  are  obtained: 


i 

I 


i 


l 

I 


i 

| 

t 


1 


x(t)  =  f  (x  (t)  +  fix  (t)  ,t) 
o 


(6) 


or. 


x(t)  =  f  (xQ(t) ,t)  +  A (t) 


5 x  ( t )  +  H.O.T. 


(7) 


x  (t) 
o 


where:  A(t) 


x  (t) 
o 


=  _3_f 
3x 


is  the  matrix  of  partial 


x  (t) 

o 


derivatives  with  respect  to  the  state  variables  evaluated  along  x  (tl , 

o 

the  reference  solution. 

Subtracting  Equation  5  from  Equation  7  on  both  sides  and 
neglecting  the  higher  order  terms  yields  the  basic  perturbation  state 
relationship  of  the  differential  corrector: 


6x(jt)  =  A  (t) 


5x(t) 


(8) 


x  (t) 
o 


Since  A(t) 


is  evaluated  along  the  reference  solution 


x  (t) 
o 


trajectory,  one  generally  obtains  a  time  dependent  linear  differential 
equation  which  has  the  solution: 


6x(t)  =  <J>(t,t  )  6x(t  ) 
o  o 


The  state  transition  matrix  4>(t,t  ),  is  obtained  by  solution  of: 

o 


(9) 


2Y 


- :  _ 


d  U(t,t  )]  =  A(t) 

dl 


v(t,t  ) 

o 


x  (t) 
o 


with  initial  conditions,  'i(t  ,t  )  =  I,  the  identity  matrix. 

o  o 

The  general  observation  relationship  of  Equation  3  yields  a 
similar  linear  perturbation  relationship  between  the  changes  in  the 
observations  and  the  state  variables.  Evaluating  the  observation 
equation  along  the  reference  trajectory,  x  (t)  ,  at  discrete  times,  t^, 
yields  the  nominal  measurements: 


z  (t  )  =  h  (x  (t  )  ,t  ) 
on  o  n  n 


which  is  defined  without  a  corruptive  noise  term.  Note  that  a  random 

noise  corruption  term  is  included  on  the  actual  observations,  z(t  ). 

n 

These  true  observations  satisfy: 


z(t  )  =  h  (x  (t  )  +  <5x (t  ),t  )  +  r  (t  ) 
n  on  n  n  n 


Expanding  about  the  reference  solution  in  a  Taylor's  series 


yields : 


z(t  )  =  h  (x  (t  ),t  )  +  H  (x  (t  ),t  )  6x(t  )  +  H.O.T.  +  r(t  )  (13) 

n  onn  onnn  n 

Subtracting  Equation  11  from  both  sides  of  Equation  13  yields  a  rela¬ 
tionship  for  the  residuals  of  the  observations,  v(t  )  -  the  difference 

n 

between  the  true  and  reference  observations: 


v ( t  )  5  z(t  )  -  z  (t  )  =  H  (x  (t  ),t  )  6x(t  )  +  r(t  ) 
n  non  onn  n  n 


where:  i)  the  H.O.T.  have  been  neglected,  and  ii)  H(x  (t  ),t  )  is  the 

on  n 

partial  derivative  matrix  of  the  geometry  relationship  evaluated  at 


the  observation  times  along  the  reference  trajectory.  Note  that  the 


observation  residuals  of  Equation  14  differ  from  what  is  often  referred 

to  as  the  estimator  residuals  -  the  difference  between  the  predicted 

and  the  estimated  state  solutions. 

Assuming  the  validity  of  the  linearization  process,  one  now  has 

a  linear  dynamics  and  geometry  relationship  from  which  the  differential 

corrector  may  be  developed.  To  initiate  this  process,  one  seeks  to 

minimize  the  difference  between  the  measured  observations,  z(t  ),  and 

n 

the  set  of  reference  observations,  z  (t  )  ,  in  a  cost  function  which  is 

o  n 

a  weighted  quadratic  form.  Expanded  to  first  order,  the  cost  function 
is : 

J  =  I i  z<cn)  -  2Q(tn)  -  H(xo(tn),tn)  «x(tn)  ||2  (15) 

W 

which  in  general  is  weighted  by  the  weighting  matrix,  W. 

As  Reference  13  shows,  if  W  is  chosen  to  be  the  positive 

definite  inverse  observation  covariance  matrix,  R(t  )“ 1 ,  the 

n 

minimization  of  J  results  in  a  first  order  approximation  to  a  minimum 

variance  estimate  for  the  state  variables.  For  this  application, 

R(t  )  is  chosen  as  diagonal  under  the  conditions  that  the  random  errors 
n 

in  azimuth  and  elevation  are  uncorrelated  within  a  given  observation. 

The  estimate  is  shown  to  be  unbiased  (13)  under  the 
restrictive  assumptions  of: 

-  Exact  dynamics 

-  Zero  mean,  random  observations 

-  Valid  linearization  assumptions 

The  numerical  results  of  Section  II. D.  will  show  the  efficacy  of  these 
assumptions  for  the  current  reentry  application. 
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In  actual  application,  the  batch  processor  seeks  to  find  the 

6x(t  )  correction  at  a  selected  epoch,  t  ,  which  will  minimize  the  cost 
o  o 

function.  For  structuring  this  batch  processor,  it  is  more  convenient 

to  formulate  the  above  problem  in  matrix-vector  form  which  incorporates 

all  n  observations  to  be  processed  by  the  estimator  for  a  given  epoch 

update.  The  notation  below  includes  capital  letters  for  matrices, 

with  subscripts  reflecting  the  observation  index  being  processed,  T  , 

n 

and  parenthetical  subscripts  for  matrices  including  a  collection  of  up 

to  n  such  observations,  T. 

W 

Define , 


T 


(n) 


and , 


r 


(n) 


H(x  (t  )  ,t  )  ' (t  ,t  ) 
o  n  n  no 


(16) 


V  _  / 

H(xo(ti) ’tf)  t(t1»to) 


f  r(tn}  \ 

^n-l> 

<  ’  >  (1?) 
'  ^(tx)  ' 


The  matrix  equivalent  to  Equation  14  then  becomes. 
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A 


satn.*- 


... 


(18) 


v.  .=  T,  .  Sx(t  )  +  x .  . 
(n)  (n)  o  (n) 


Define  R.  .  as  the  matrix  of  the  observation  noise  covariances  in  a 

(n) 


symmetric,  block  diagonal  form: 


(n) 


f  R(tn> 


R(tn-1> 


(19) 


R(t1)  ) 


R^  is  block  diagonal  under  the  conditions  that  the  random  error  in 

successive  observations  (ie,  at  times  t  ,  and  t  )  is  uncorrelated. 

n-i  n 

This  is  a  reasonable  assumption  for  the  10  second  interval  between 
observations. 

The  weighted  least-squares  cost  function  J  is  now  of  the  form: 


J  ‘  [v(»)  -  T(n)  %)  1*(„)  -  T«0  Sx(to)]  (20) 

Since  is  positive  definite,  the  necessary  and  sufficient 


condition  for  minimizing  J  is: 


3  J 


=  0T 


36x(t  ) 
o 


(21) 


Solving  for  6x(t  )  yields  (13): 
o 


“V-'Vl'V'V"’*.)1  R(n)";(n) 


(22) 


By  definition,  for  zero  mean  6x(t  ),  the  state  covariance  matrix  for 

o 

the  differential  corrector  is: 
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S  (t  )  =  E  (6x(t  )  6x(t  )  ) 
n  ,n  o  o  o 


S  (t  )  denotes  the  state  covariance  matrix,  and  E  is  the  expectation 
n,n  o  r 

operator.  Substituting  <5x (t  )  f rom  Equation  22  yields: 

Sn,n(to>  =E  ^<T(n)T  R(n)“  1  V'  "(n/Vf1  %)  * 

^n))"1!  <24) 

"  (T(n)TR(n)'lT(n)>‘lT(n)TR(n)-lE^(nr(n)T>  * 


Rr  -1  T  .  .  (T  .  ,T  R .  ”  1  T.  .)' 
(n)  (n)  (n)  (n)  (n) 


Recall  that  v,  ,  is  the  accumulated  vector  of  the  observation 
(n) 

residuals.  By  definition,  for  zero  mean  observations: 

%)  S  «;W;(n)T>  (26> 

Therefore,  after  substitution  and  simplification,  the  state  covariance 
expression  becomes: 


S„,„(to>  •  «(„)T  *(„)“  T(„>>“ 


Should  a  priori  state  covariance  information  (denoted  here  by  S'  (tQ)) 
also  be  available,  it  can  be  shown  that  (14): 


S  (t  )  =  (S'(t  )_1  +  T,  .  R.  ~1  T,  ,)-1 
n,n  o  o  (n)  (n)  (n) 


Assuming  a  constant  ballistic  parameter,  3,  over  some  portion  of 
the  trajectory,  the  application  of  the  differential  corrector  reduces 
to  the  processing  of  a  series  of  segments  of  observations.  Each 
update  epoch  may  be  designated  with  the  index  m.  At  an  epoch,  the 


X  I 


estimator  may  process  up  to  k  observations  for  an  update.  In  the 
limiting  case,  only  one  observation  per  epoch  update  will  be 
processed.  The  figure  below  illustrates  the  notation  for  the 
sequential  estimator  formulation.  The  state  estimate  at  epoch  time 


t  is  denoted  x  ,(t  )  prior  to  update.  At  time  t  ,  the  estimator 
m  m-1  m  m 


may  process  one  or  more  observations.  The  updated  state  estimate. 


xm(t  ) ,  is  obtained  and  is  then  propagated  forward  to  the  new  epoch 


time,  t  ^ .  The  process  then  repeats. 


m-l^m-* 


t  n=2 
m 

n 


t  n=k 

•  m 

w  n 


Observations 


t  ..  n=2 


) Observations 


Having  propagated  forward  from  an  epoch  at  time  t  ,  the  a  priori 

state  covariance  at  the  new  epoch,  t  is  S  .  Using  the 

m+1  m+1 ,m 

reference  trajectory,  xm(t)  propagated  forward  by  integration  of  the 
dynamics  differential  equation,  one  formulates  the  information  and 
weighted  residuals  matrices: 


(T.  ,T  R,  ~  1  T,  ,)  and  T.  ,T  R.  "  1  v ,  , 
(n)  (n)  (nr  (n)  (n)  (n) 


where  n  =  1,  .  .  .  ,  k,  the  number  of  observations  used  to  update  the 


state  estimate  at  the  epoch  t  . 

m 

An  initial  correction  to  the  state  epoch  is  then  found  from: 


6x(tm+l)  = 


(S 


-  i 


m+1  ,m 


+  T 


(n)  *(n) 


-  i 


T(n)> 


-  x 


(n)  R(n) 


-  l 


(n) 


(29) 


This  6x(tm+.;  is  added  to  the  reference  trajectory  to  form  a  new 
reference  trajectory  in  an  iterative  fashion,  such  as: 


m+1 


- 


<tm+l> 


6x(tm+l>i 


(30) 


Successive  application  of  the  differential  corrector  continues 
until  the  convergence  criteria  on  the  observation  residuals  are 
satisfied.  The  initial,  a  priori  state  covariance  matrix  remains 
constant  until  the  iterative  process  has  converged.  On  the  final 
iteration,  the  following  update  results: 


Xm+1  (tm+l)  ”  Xm  ^m+l*  +  ^  5x(tm+l)i 


Sm+1  ,m+l  "  (Sm+l,m_1  +  T(n)T  R(n)_1  T(n)0) 


(31) 

(32) 


where  l  equals  the  number  of  iterations  required  for  convergence.  The 


information  matrix,  T 


(n)  (n)  (n) 


T  .  .  ,  is  evaluated  from  the  reference 


trajectory  on  the  last  iteration. 

For  orbital  applications,  the  standard  convergence  criterion  of 

the  iterative  process  has  been  implemented  by  a  magnitude  check  on  the 

observation  residuals  (5,13).  One  may  define  a  norm  of  the  observation 

residuals  within  \  .  .  on  a  given  iteration  as: 

(n) 

k  p  i 

ilv.ll  5  [  Z  (  I  v  2)  ]*  (33) 

i-1  j=l  2 

indexing  over  j  for  the  components  within  a  single  observation  residual 
vector,  and  over  i  for  all  the  observations  being  processed  for  the 
current  epoch  update. 

The  convergence  criteria  have  often  been  specified  in  terms  of 
either  a  relative  or  absolute  criterion  (5,13): 


'VI  -  '  ’  vi+l ' 


<  ci  Relative 


jv^l |  <  e 2  Absolute 


When  either  of  these  conditions  are  satisfied,  the  iterative  update  is 


considered  to  have  converged.  Note  that  Equations  34  and  35  do  not 
represent  a  sufficiency  condition  for  convergence  of  an  iterative 


process.  Kaper,  et  al  (15)  point  out  that  the  above  criterion  only 
satisfy  a  necessary  condition  for  convergence  of  an  iterative  process 


of  successive  linearizations  for  a  nonlinear  system.  They  suggest 
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that  a  "sufficiency  criterion"  must  also  be  specified  which  examines  the 
magnitude  of  the  resulting  state  correction,  on  the  final  iteration. 

This  <5x.(t  i  =  1,...,8  must  also  be  small  for  the  process  to  have 

converged.  In  actual  estimator  application,  the  e*  and  C2  necessary 
conditions  were  chosen  empirically  such  that  the  state  correction  on 
the  converged  iteration  yielded  values  less  than  approximately  .01 
times  their  respective  standard  deviations  from  the  estimator-computed 
state  covariance  matrix.  While  not  a  true  sufficiency  condition  in  a 
strict  mathematical  sense,  this  allows  implementation  on  the  digital 
computer  without  generation  of  many  additional  iterations  of  the 
estimator,  whose  solutions  contribute  little  to  modifying  the  state 
estimate. 

To  summarize  the  conditions  used  to  determine  convergence: 

1)  Assure  that  the  necessary  criterion  of  the  observation 
residual  magnitude  of  either  Equations  34  or  35  is  satisfied. 

2)  Assure  that  the  individual  6x^  on  the  final  iteration,  l , 

are  less  than  approximately  .01  times  the  standard  deviation 

of  their  respective  terms. 

i 

Additionally,  one  should  examine  the  following  conditions: 

3)  Assure  that  the  absolute  value  of  the  residuals  is  less  than 
the  one  sigma  magnitude  of  the  observation  covariance: 

|v  |  <  /~R  i=l  azimuth  i=2  elevation  (36) 

1  _  nii 

Failure  of  the  bound  of  Equation  36  indicates  a  potential 
onset  of  divergence. 

4)  Assure  that  the  total  change  in  each  state  variable  is  less 
than  or  equal  to  its  one  sigma  a  priori  uncertainty: 


£  6x .  .  <  /  S  “ 

.  .  lj  -  nrt-1  ,m.  . 

J=1  11 


i  =  1,8 


(37) 


This  last  criterion  provides  assurance  that  the  state  upd^.e  at 
a  given  epoch  is  consistent  with  the  state  uncertainties 
defined  by  the  use  of  ^  1  within  the  update  expr  'sion  for 

6x(t^_^)  as  the  dynamics  weighting  matrix. 

As  later  discussions  indicate,  violation  of  this  criterion  offers 
a  valuable  measure  of  the  potential  onset  of  divergence  of  the  state 
solution,  for  the  non-exact  dynamics  cases.  In  the  circumstance  where 
observation  residuals  do  not  grow  unbounded,  or  are  not  non-zero  mean, 
this  criterion  may  offer  a  measure  which  indicates  a  need  to  implement 
a  fading  memory  or  other  compensation  method  to  the  estimatcr. 

Application  of  this  batch  processing  algorithm  to  reentry  esti¬ 
mations  has  often  resulted  in  poor  estimator  performance.  This  is 
largely  due  to  the  more  significant  non-linear  dynamics  of  reentry  and 
the  use  of  a  deterministic  and  simplistic  dynamics  model  of  Equation  2 
in  an  uncertain  dynamics  region.  Consider  the  effects  of  neglected 
nonlinearities.  Both  dynamic  and  geometric  nonlinearities  can  impact 
the  state  vector  and  covariance  update.  As  a  means  to  insure  the 
validity  of  the  linearization  assumptions,  a  linearity  check  was 
utilized  which  was  applied  to  each  iteration  of  the  differential 
corrector  process  (16).  This  was  done  to  insure  reasonable  confidence 
in  the  linear  relationship  between  changes  in  the  state  variables  and 
the  observations,  as  shown  in  Equation  14. 

Consider  the  reference  trajectory  observations  on  two  successive 


3L 


-  - 


iterations,  i  and  i+1: 


(38) 


2  (t  ) .  =  h(x  (t  ) . ,t  ) 
n  m  1  m  m  l  m 

n  n  n 


2  )i4.,  =  h(x  (t  )  ,  1 1  *  b  ) 

n  m  i-rl  m  m  I+l  m 

n  n  n 


(39) 


Assuming, 


x  (f),  =  X  (t)  +  <$x(t) 

m  i+l  m  l 


(40) 


and  expanding  about  the  ith  reference  trajectory  yields: 


z  (t  ).,,.  =  h(x  (t  ).,t  )  +  H(x  (t  )  ,t  )  5x(t  )  +  H.O.T.  (41) 

n  m  i+l  mmim  tnmim  m 

n  n  n  n  n  n 


Therefore • 


6z  =  z  (t  )  ,  -  z  (t  )  =  H(x  (t  ),,t  )  <5x(t  )  +  H.O.T.  (42) 

n  n  m  i+i  nmi  mmim  m 

n  n  n  n  n 

or  as  implemented  relative  to  the  epoch  update  of  the  state  variables: 


6z  =  H(x  (t  )..t  )<J>(t  ,t  )  6x(t  )  +  H.O.T. 

n  mmim  mm  m 

n  n  n 


(43) 


One  may  specify  a  H.O.T.  magnitude  check  as: 


H.O.T. 


lT(n)^<^ 


<  £ 


(44) 


on  each  successive  iteration  of  the  iterative  solution  process.  The 

H.O.T.  is  obtained  from  Equation  43,  since  5z  and 

n 

H(x  (t  ).,t  )  <p(t  ,t  )  are  known  quantities, 

mmim  mm 

n  n  n 

In  actual  application,  e  was  chosen  as  0.1  since  this  lineariza¬ 
tion  retains  the  first  term  in  the  Taylor's  series  expansion.  This 
linear  term  should  be  at  least  one  order  of  magnitude  greater  than  the 
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H.O.T.  in  the  expansion.  If  this  criterion  is  violated,  one  could 


abandon  the  linearized  approach  to  estimation  with  a  deterministic 
dynamics  model  or  reduce  the  time  span  over  which  the  linearization  is 
assumed  valid.  Limiting  the  number  of  observations  processed  at  a 
given  epoch  update  limits  the  time  span  over  which  the  linearization 
approximation  is  valid. 

With  exact  dynamics,  the  later  Monte  Carlo  results  (Figures  6  - 
13)  illustrate  the  basic  validity  of  the  estimator  performance  using 
this  Taylor's  series  linearization  approach.  As  the  non-exact  dynamics 
cases  substantiate,  the  dominant  impact  on  estimator  performance  is 
due  to  the  uncertainties  in  the  true  reentry  dynamics.  This  circum¬ 
stance  causes  one  to  address  the  dynamic  uncertainties  first,  prior  to 
pursuit  of  more  sophisticated  techniques  such  as  a  higher  order 
estimator  or  a  statistical  linearization  approach.  These  points  are 
discussed  in  Chapter  III. 

To  summarize,  the  estimator  update  and  propagation  equations  are 
provided  below  for  an  infinite  memory  formulation: 

Propagations  between  epochs: 

State:  integrate,  with  initial  conditions  x  (t  ),  from  the 

m  m 

epoch  at  t 

m 

x(t)  =  f (x(t) ,t)  (45) 

to  obtain  x  (t  ,,)  at  the  epoch  t 

m  m+1  r  m+1 

Covariance: 


m+1  ,m 


S„  „  ♦<»_,,  .O' 

m,m  m+i.  m 


(46) 
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Update  at  next  epoch: 
State : 


WW  =  VW  +  Z  6x(tmfl)i 

i=l 


(47) 


Covariance : 


C  =  /C  —  1  4.  T  13  -  1  T  1 

ra+l,m+l  m+1  ,m  (n)  (n)  (n)^ 


(48) 


Recall  that  SL  is  the  number  of  iterations  required  to  satisfy  conver- 
T 

gence  and  T,  N  R,  1  T,  .  is  evaluated  from  the  reference  trajectory 
(n)  (n)  (n)  J  J 

on  the  final  iteration,  l. 
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B.  Dynamics  Model 

B.l.  Eight  Dimensional  Formulation 

Previously  stated,  an  exponential  atmosphere  and  constant 
ballistic  parameter  assumption  were  combined  with  the  position  and 
velocity  variables  to  construct  an  eight  dimensional  state  vector, 
x(t).  In  this  manner,  the  estimator  solution  will  contain  two 
distinctly  quantifiable,  physical  parameters:  the  ballistic 
parameter,  6,  and  the  atmospheric  density  scale  height,  Q. 

Consider  the  influence  of  the  atmospheric  density  model  on  the 
estimator  performance.  A  review  of  the  standard  atmosphere  density 
equations  illustrates  an  increased  complexity  when  compared  to  an 
exponential  (or  isothermal)  density  model.  The  application  of  the 
standard  density  model  requires  separate  equations  in  several 
altitude  layers  during  reentry.  All  densities  are  defined  relative  to 
base  altitude  values  for  density,  p  ,  and  temperature,  T  .  In  the  non- 
isothermal  layers  the  density  is  functionally  related  to  the  direction 
and  rate  of  change  of  the  temperature  with  altitude.  This  relationship 
is  characterized  by  the  thermal  lapse  rate,  L.  Equation  49  defines  the 
non-isothermal  standard  density  expression. 


Tb  +  L(H*  -  H*b) 


{1  +  (Mg  */RL)  } 
o  o 


(49) 


where : 


H  ,H  =  geodetic  altitude,  and  geodetic  base  altitude, 
B 

respect ively 

M  =  molecular  weight  of  air 

o 


J5 


g  *  =  acceleration  of  gravity 

o 

R  =  universal  gas  constant 

M  and  R  are  generally  constant  below  80  KM  altitude  where  continuum  gas 
dynamics  apply,  although  the  standard  atmosphere  equations  are  often 
used  well  above  86  KM  altitudes.  (4,10) 

Equation  50  defines  the  isothermal  standard  density  expression. 

Mg*  (H*  -  H*  ) 
o  o  B 

p  =  p  exp  { - }  (50) 

B 

Within  the  linearized  estimation  technique  one  must  linearize  the 

equations  of  motion  relative  to  a  reference  trajectory.  This  involves  a 

Taylor’s  series  expansion  about  the  reference  trajectory  solution.  This 

requires  partial  derivatives  of  the  density  expressions  with  respect  to 

the  trajectory  position  variables  (x,y,z).  The  position  variables  are 

* 

embodied  within  the  geopotential  altitude  terms,  H  .  Due  to  the 
* 

presence  of  H  within  the  denominator  of  Equation  49,  the  complexity 

and  nonlinear  character  of  the  partial  derivatives  are  more  pronounced 

within  the  non-isothermal  layers  of  the  atmosphere. 

The  derivatives  of  the  density  expression  with  respect  to  x,  y 

and  z  are  also  difficult  to  calculate  across  the  altitude  layers  of  a 

standard  density  model.  The  linearization  is  often  required  to  be 

valid  across  layers  of  differing  thermal  lapse  rate  or  between  non- 

isothermal  and  isothermal  layers.  This  will  impact  the  linearization 

process  since  the  density  across  these  layer  boundaries  is  a 

* 

continuous  function  of  H  ,  but  the  spatial  first  derivatives  of  the 
density  are  not  continuous.  The  alterate  formulation  of  an  isothermal 
density  model  for  the  entire  atmosphere  can  greatly  simplify  the 
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implementation  aspects  of  the  estimation  technique. 

* 

P  =  PQ  exp  {-  H  /Q}  (51) 

where : 

p  =  Earth  surface  reference  density 
o  J 

•k 

H  =  geodetic  altitude 

Q  =  RT  is  defined  as  a  constant  in  the  altitudes  of 

o 

Mg* 

o  o  reentry,  called  the  density  scale  height. 

* 

(This  is  equivalent  to  Equation  50  if  T  =T  ,  H  =0,  and  p  =p  .) 

B  o  B  Bo 

While  this  exponential  model  is  an  inherently  poorer  functional 
description  of  the  "mean"  atmosphere  density  variations  with  altitude, 
it  is  more  easily  incorporated  into  the  reentry  estimator.  If 
acceptable  estimator  performance  can  be  available  from  this  more  sim¬ 
plified  expression,  it  will  greatly  reduce  the  complexities  associated 
with  use  of  the  standard  atmosphere  equations.  The  advantage  of  this 
simplified  implementation  is  a  reduced  mathematical  complexity  and  the 
availability  of  continuous  valued  density  and  partial  derivatives  of 
density  along  the  reentry  trajectory. 

The  density  scale  height  was  chosen  as  an  eighth  state  variable 
since  it  is  slowly  changing  within  the  altitudes  of  reentry  (below  100 
KM)  as  shown  in  Figure  2.  It  also  is  the  basic  independent  variable 
within  the  exponential  atmosphere  density  expression  (Equation  51). 

Alternate  state  variable  forms  were  considered  to  augment  the 
basic  six  position  and  velocity  variables.  An  approach  to  estimating 
the  uncertain  dynamic  variables  has  been  used  by  many  authors  which 
incorporates  a  first  order  Gauss-Markov  assumption  for  the  unmodeled 
accelerations.  Tapley  and  Ingram  (17)  suggest  such  a  first  order 
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Gauss-Markov  process.  In  their  formulation,  the  following  expression 


is  used: 

e^(t)  =  -B^  e  (t)  ^  +  u^(t)  i  =  1,2,3  (52) 

where : 

e^  is  a  time  correlated  component 
are  constant  coefficients 

and,  u^(t)  is  a  white,  Gaussian  noise  term  with  statistics: 

E  [u(t) ]  =  0  (53) 

and , 

E  [u(t),u(t  +  t)T]  =  Q (t)  6 (t)  (54) 

Similarly,  Myers  and  Tapley  (18)  formulate  the  unmodeled  accelerations 
as : 

e±(t)  =  -  e^O/T.  +  u.(t)  (55) 

where  T.  =  a  time  correlation  coefficient  =  1/B..  This  approach  con- 

l  l 

tains  several  drawbacks  for  the  general  decay  satellite  application. 
First,  it  requires  a  pseudo-noise  compensation  to  the  model  dynamics, 
while  the  deterministic  estimator  dynamics  of  Equation  2  do  not.  The 
time  correlation  coefficients  must  be  determined  along  with  the  Q(t) 
covariance  matrix  elements,  requiring  extensive  simulation  analysis 
for  tuning.  As  the  orbital  experience  indicates  (5)  ,  this  process  does 
not  have  high  potential  for  success  in  the  uncertain  reentry  dynamics 
regions  with  short  arcs  of  empirical  data  to  assist  the  tuning  process. 
Secondly,  if  separately  estimated  along  with  the  e^  terms,  the  T^  or  B^ 
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coefficients  can  ultimately  increase  the  state  dimension  from  the 
original  six  positions  and  velocities  to  12  state  variables.  Lastly, 
as  Myers  and  Tapley  (18)  point  out,  the  specification  of  initial 
conditions  and  covariance  values  for  these  variables  is  not  easily 
related  to  the  physical  processes  present  in  reentry.  Minor  errors 
in  specification  of  these  initial  values  were  found  to  have  signifi¬ 
cant  impact  on  the  estimator  performance. 

Initial  values  and  covariance  magnitudes  for  the  ballistic 
parameter,  8,  and  the  density  scale  height,  Q,  are  readily  available. 
The  value  of  8  and  its  associated  covariance  is  generally  available 
from  propagating  the  last  SPACETRACK  orbital  estimation  to  the  first 
reentry  epoch.  The  initial  value  of  Q  and  its  covariance  can  be 
developed  from  a  standard  atmosphere  model  and  modified  with  knowledge 
of  the  local  atmosphere  properties  which  may  possibly  be  known  for  the 
specific  reentry  case  of  interest. 

In  the  current  application,  the  initial  density  scale  height 
was  derived  from  a  least  squares  fit  to  the  base  density  values  of  the 
altitude  layers  within  the  U.S.  1962  Standard  Atmosphere  (Appendix 
B.I.).  The  resulting  scale  height  was  Q  =  7.0031  KM. 

A  rectangular,  inertial  coordinate  frame  (Earth  Centered 
Inertial  -  ECI)  was  chosen  to  minimize  the  complexities  of  the 
computational  procedure  (19).  As  such,  the  following  state  vector  was 


selected : 


X 


(56) 


The  product  Bp  was  chosen  so  that  x_  would  be  of  order  one  in  the 
o  7 

kilometer,  kilogram,  second  (KKS)  system  used  for  the  dynamic  computa¬ 
tions.  Use  of  the  more  standard  meter,  kilogram,  second  (MKS)  system 
would  result  in  very  small  covariance  terms  for  the  x,  state  (Bp  )  and 
aggravate  an  already  ill-conditioned  state  covariance  matrix  (to  be 

shown  later  in  Section  D.3.).  The  term  p  equals  the  Earth  surface 

o 

reference  density  model  value,  and  Q  is  the  density  scale  height  of  the 
exponential  density  expression: 

★ 


P 


e”  (H  /Q> 


(57) 


The  geodetic  altitude,  H  ,  is  developed  from  the  Earth  model  parameters 
as  follows: 


H*  =  g/g0Ho  (58) 

where:  g  is  local  gravitational  acceleration  obtained  from  the 

magnitude  of  the  geopotential  acceleration  components  along 
the  x,  y,  and  z  coordinates. 

gQ  is  the  reference  geoid  level  gravitational  acceleration 


.  A.  y 
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The  Ho  term  represents  a  geocentric  altitude: 

H  =  R  -  R  (59) 

o  s 

R  is  the  local  radius  position  relative  to  the  Earth  center, 

R  =  (x2  +  y2  +  z2)^,  and  Rg  is  the  local  Earth  surface  radius  position: 

R  =  R  (1  -  f)  [1  -  (2f  -  f2)  cos2  6]  **  (60) 

s  o 

where : 

f  =  the  flattening  factor  of  the  reference  geoid  whose 

elliptical  shape  is  consistent  with  the  dominant  non¬ 
central  gravity  term,  J^.  due  to  the  equitorial  bulge. 

6  =  local  latitude,  cos  6  =  (x2  +  y2)^/RQ 

R  =  radius  of  the  reference  geoid  at  the  equator 
o 

The  estimator  dynamics  model  includes  both  aerodynamic  and 
geopotential  acceleration  terms  within  the  f(x(t),t)  expression.  In 
this  manner,  the  dynamics  equation  will  have  the  form: 


(61) 


The  aerodynamic  accelerations  are  derived  from: 
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ad  =  ^  0P  VR  \ 


where  V  is  Che  scalar  magnitude  of  the  velocity  vector  relative  to  a 
rotating  atmosphere.  The  decomposition  of  this  aerodynamic  accelera¬ 
tion  along  the  three  coordinates  x,  y  and  z  results  in: 

f d  =  -%■  &P  VR  (x  +  wy) 

x 

fd  =  ~h  Bp  VR  (y  -  'ox)  (63) 

y 

fd  =  ^  6p  VR  * 
z 

where  w  is  the  scalar  magnitude  of  the  rotational  velocity  of  the  Earth, 
about  the  z  axis. 

The  following  development  yields  the  geopotential  accelerations 
(20).  Consider  the  geopotential  function: 


n  GM  R 

CO  , 

=  E  l  - r 


n=0  m=0  R 


~~ —  P®  (6')  [C  cos  (mX)  +  S  sin  (mX)J 
n+1  n  nm  nm 


where:  G  =  universal  gravitational  constant 

M  =  mass  of  Earth 

Rq  =  mean  equitorial  Earth  radius 

R  =  distance  from  Earth  center 

=  associated  Legendre  functions  of  degree  n  and  order  ra 

C  ,S  =  zonal,  tesseral  harmonic  coefficients 
nm  nm 

6'  =  z/R 
X  =  longitude 

One  may  define  a  zonal  and  tesseral  harmonic  component  such  that: 


um  =  — -rr5-  Pm  (<5 ' )  cos  (mX) 
n  n+1  n 

K 


A3 


GM  R 

V™  =  - T\ -  Pm  (61)  sin  (mX) 

n  n+1  n 

K 


(66) 


The  geopotential  which  results  is: 


n 

00  mm 

<J>  =  E  E  (C  U  +  S  V  ) 
^  run  n  nm  n 
n=0  m=0 


(67) 


Defining  ag  as  the  Right  Ascension  (RA)  of  the  Greenwich 
Meridian,  one  obtains  the  following  geopotential  acceleration  terms: 


f  =  cos  a.  VV  sin  a  VV 
g„  s  x  s  y 


f  «  sin  a  VV  +  cos  a  VV 
gy  xx  s  y 


(68) 


f  =  VV 
gz 


where  the  gradient  terms  are  defined  by: 


CO 

=  E 
n =0 

n 

E 

m=0 

(c 

nm 

aum 

n 

s 

nm 

av”1 

n 

3x 

T 

3x 

oo 

=  E 
n=0 

n 
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(69) 


The  Smithsonian  Astrophysical  Observatory  SAO-III  (7)  Earth 
Model  was  used  to  define  the  model  parameters  and  zonal  and  tesseral 
coefficients.  In  actual  estimator  development  and  simulation  only  the 
central  gravity  Cqq  term  and  Earth  oblateness  ^20^2^  term  were  used  in 
the  reentry  altitudes  regions.  The  computer  program  code  is  easily 
modified  by  changing  dimensions  and  indexing  to  incorporate  the  full 
SAO-III  Earth  Model.  The  retention  of  the  component,  and  no  higher 
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order  geopotential  terms,  is  consistent  with  the  relative  geopotential 


and  aerodynamic  acceleration  magnitudes  and  uncertainties  in  the  reentry 
altitude  regions.  An  ability  still  exists  to  use  the  full  SAO-III 
Earth  Model  in  propagating  the  final  orbital  vector  to  the  first  epoch 
update  point  of  the  reentry  observations.  Consider  the  relative  aero¬ 
dynamic  and  geopotential  acceleration  magnitudes. 

1)  In  a  high  altitude  region  above  the  reentry  data 

points  (approximately  60  NM  or  111  KM  high) : 

Aerodynamic  acceleration:  a^  £  5.75  x  10“ 13  K/s2 

Central  gravity  term:  a^  =  .00947  K/s2 

o 

J0  gravity  term:  aT  =  1.48  x  10” 5  K/s2 
J2 

J.  gravity  term:  a.  =  4.32  x  10  8  K/s2 
4  J4 

This  illustrates  the  need  to  carry  the  higher  order 

geopotential  terms  when  propagating  the  final  orbital 

solution  forward  to  the  first  reentry  data  point.  In 

this  region,  the  higher  order  geopotential  terms  yield 

accelerations  whose  magnitudes  are  significant  relative 

to  the  magnitude  of  the  aerodynamic  acceleration.  (All 

harmonics  of  order  or  higher  exhibit  accelerations 

whose  magnitudes  are  similar  to  or  smaller  than  the 

J,  term.) 

4 

2)  In  the  reentry  altitude  region  (IR  data) : 

(appioximately  40  NM  or  74  KM  high) 

Aerodynamic  acceleration:  £  5.41  x  10”5  K/s2 

Central  gravity  term:  aj  =  .00957  K/s2 

o 

gravity  term:  a.  s  1.52  x  10-5  K/s2 
L  J2 


i 


J,  gravity  term:  a  =  7.4  x  10-8  K/s2 
J4 

The  J2  acceleration  contribution  becomes  of  the  same  order 
of  magnitude  as  the  aerodynamic  acceleration.  It  is  there¬ 
fore  necessary  to  retain  the  term  in  the  dynamics  model. 
The  higher  order  terms  may  be  neglected. 

3)  Typical  Earth  impact  accelerations: 

(assume  Vs  .15  K/s) 

K 

Aerodynamic  acceleration:  a^  s  5.5  x  10-3  K/s2 

Central  gravity  term:  a^  =  .009798  K/s2 

o 

J„  gravity  term:  a  =  1.59  x  10~5  K/s2 
L  J2 

The  aerodynamic  accelerations  dominate  the  higher  order 
geopotential  terms  and  assume  magnitudes  of  an  order  similar 
to  the  central  gravity  accelerations.  The  J 2  gravity  term 
accelerations  become  negligible  relative  to  the  central 
gravity  term  and  aerodynamic  accelerations. 

A  similar  comparison  can  be  made  between  the  acceleration 
magnitudes  of  the  individual  geopotential  terms  and  the  anticipated 
uncertainties  within  the  aerodynamic  accelerations  of  a  simulated  true 
reentry  dynamics  set.  Using  an  example  from  a  typical  single  run  of 
the  Monte  Carlo  analysis  (with  a  truth  model  containing  the  U.S.  1962 
Standard  Atmosphere  and  8  as  a  function  of  Mach  number)  yields  the 
following  two  cases: 

CASE  1:  Assume  x,  y,  z,  Sp^,  Q  to  be  one  standard  deviation 
from  their  mean  estimator  determined  values 
CASE  2:  Assume  only  x,  y,  z  to  differ  by  one  standard  deviation 
from  their  mean  estimator  determined  values 


The  deviations  from  the  mean  aerodynamic  accelerations  at  a  high 


and  low  reentry  altitude  region  of  IR  data  are: 

73.04  KM  altitude  - 

CASE  1:  AaD  =  3.66  x  10"5  K/s2 

CASE  2:  =  5.466  x  10"6  K/s2 

41.23  KM  altitude  - 

CASE  1:  Aa^  =  2.346  x  10"3  K/s2 

CASE  2:  Aa^  =  2.748  x  10'*4  K/s2 

In  the  beginning  of  the  IR  data  region  (73.04  KM),  the  aero¬ 
dynamic  uncertainties  can  be  of  the  same  order  as  the  acceleration 
magnitude.  During  the  reentry  data  region  (41.23  KM),  the  uncertain¬ 
ties  are  of  the  same  order  or  much  greater  than  the  acceleration 
magnitude  depending  on  the  random  error  in  each  of  the  pertinent  state 
variables.  They  can  even  become  significant  relative  to  the  magnitude 
of  the  central  gravity  accelerations  (approximately  .0096  -  .0098  K/s2). 
As  a  result  of  these  relative  acceleration  considerations,  the  estima¬ 
tor  model  used  both  the  central  gravity  and  Earth  oblateness  (J2)  terms 
from  the  geopotential  expansion. 

B . 2.  Seven  Dimensional  Formulation 

As  previously  mentioned,  the  standard  practice  by  operational 
agencies  such  as  the  USAF  SPACETRACK  System  (5)  has  been  to  use  a  seven 
dimensional  state  vector  for  estimator  formulation.  The  atmospheric 
density  ratios  are  ther.  obtained  from  the  standard  density  model  of 
Equations  49  and  50.  The  state  vector  becomes: 


—  - _ 
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(70) 


The  U.S.  1962  Standard  Atmosphere  Model  (10)  was  used  for  this 
density  information  with  the  geopotential  altitude  values  modified  to 
reflect  the  SAO=III  Earth  model  values  (see  Appendix  B.l).  Following 
the  same  development  as  the  eight  dimensional  estimator  formulation,  a 
complete  seven  dimensional  estimator  was  constructed.  Performance 
implications  of  each  formulation  were  made  and  will  be  discussed  with 
the  numerical  examples  in  Section  D.3.,  later  in  this  chapter. 
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C .  Observation  Relationship 

The  observation  relationship  for  the  IR  data  from  the  orbital 
sensor  takes  the  form  (see  Figure  1) : 


z(t  )  =  h(x(t  ) , t  )  = 
m  m  ra 

n  n  n 


Vt  > 

1  m 

n 


h.<t  > 

2  m 

n 


where: 


h  (t  )  =  sin“ 


(y '  z 


+  it  measured  within  the 
x'y'  plane  from  the 
negative  z'  axis  as 
in  Figure  1. 


h2(tmn)  =  Sin  1  -  (■xt^+ynrrr7TT)!s 


(measured  fr^m  the  y'z'  plane  to  the  reentry  satellite 


position  vector,  r;  see  Figure  1.) 


The  relationships  for  the  geometric  partials  matrix,  R(t  )  = 

m 

I  " 

3h(x , t) 


and  also  the  dynamics  partials  matrix,  A(t)  = 


x  (t  ) 
m  m 

n 


9f  (x,t) 


for  the  basic  eight  dimensional  estimator  are  shown 


x  (t) 
m 


in  Appendix  A. 
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D .  Numerical  Simulation  Results 

Having  developed  che  basic  estimator  structure  and  dynamics 
model,  one  must  now  quantify  its  performance  with  simulated  reentry 
data.  As  Reference  13  details,  the  application  of  such  a  successive 
linearization  technique  results  in  an  unbiased  estimate  only  under  the 
following  restrictive  assumptions: 

1)  The  inverse  covariance  matrix  R,  .-1  is  used  as  the 

(n) 

weighting  matrix  within  the  basic  cost  function  of  the 
differential  corrector  (Equation  20). 

2)  The  observations  are  randomly  distributed  and  zero  mean. 

3)  The  changes  in  the  observations  are  linear  functions  of 
the  changes  in  the  state  variables. 

4)  The  dynamics  model  is  exact. 

In  the  simulated  estimator  runs  which  follow,  the  inverse 
observation  covariance  matrix  was  used  for  weighting.  The  simulated 
observation  noise  was  generated  from  a  computer  based  random  number 
generator.  This  noise  was  essentially  zero-mean,  uncorrelated,  and 
had  a  standard  deviation  equal  to  the  standard  deviation  of  the 
observation  covariance. 

The  simulated  results  show  that  Equation  44  can  successfully  be 
used  to  limit  the  time  spans  of  updates  such  that  the  relationship 
between  observations  and  state  variables  changes  remains  essentially 
linear.  Recall  that  c  was  chosen  as  0.1  to  insure  that  the  linear 
term  of  the  Taylor's  series  expansion  was  at  least  one  order  of 
magnitude  greater  than  the  neglected  higher  order  terms. 


Particular  attention  was  given  to  non-exact  dynamics  and 
variations  in  observation  geometry  as  they  impact  the  bias  magnitude 


in  the  estimator  solutions.  The  results  of  Section  D.3.  will  also  show 


that  a  limited  observability  of  the  system  is  provided  with  the  angle 
only  observations  from  a  single  observer.  This  results  in  a  very  ill- 
conditioned  information  matrix  in  the  absence  of  a  priori  state 
covariance  data.  An  extensive  series  of  simulated  data  runs  was 
completed  with  the  basic,  infinite  memory  estimator  formulation 
(Equations  45-48).  These  were  selected  to  examine  which  factors 
required  further  evaluation  for  improvements  in  the  application. 

A  series  of  single  sample  simulated  data  runs  were  completed  to 
provide  insight  into  the  performance  of  the  estimator  while  considering 

1)  Estimator  dynamics  model  limitations 

2)  Variations  in  observation  noise  levels 

3)  Variations  in  observation  geometry  relative  to  the 
reentry  trajectory 

4)  Observations  from  multiple  observation  sources  on  the  same 
reentry  trajectory 

A  Monte  Carlo  analysis  was  then  completed  to  examine  estimator 
performance  with  exact  dynamics  at  two  different  levels  of  observation 
noise.  A  third  set  of  monte  Carlo  runs  was  completed  with  a  signifi¬ 
cant  mismatch  between  the  estimator  dynamics  and  the  dynamics  of  the 
truth  model,  from  which  the  simulated  observations  were  derived. 

The  mismatch  between  the  deterministic  estimator  and  the  truth  model 
dynamics  model  is  shown  to  be  significant.  Chapter  III  presents  a 
discussion  on  the  model  compensation  techniques  considered  to  address 
this  problem. 

Lastly,  a  number  of  special  numerical  investigations  is  pre¬ 
sented.  These  include:  a  comparison  of  the  seven  and  eight 
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dimensional  estimator  formulations,  single  observation  satellite 


observability  considerations,  and  an  analysis  of  the  propagation  of 
the  final  epoch  state  covariance  matrix  to  Earth  impact. 

D .  1 .  Single  Sample  Results 

An  extensive  series  of  single  sample  runs  were  made  to  identify 

the  most  significant  performance  aspects  of  the  basic  differential 

corrector  formulation.  In  the  application  of  the  estimator  to 

simulated  data,  the  solutions  for  x_  =  Bp  and  xQ  =  Q  are  obtained 

under  the  assumption  that  these  quantities  are  locally  constant  over 

the  time  span  between  the  epoch  update  point,  t  ,  and  the  observations 

m 

at  times  t 

m 

n 

Single  sample  simulation  results  will  be  shown  with  observations 
from  both  single  and  dual  observers  being  processed  for  each  epoch 
update.  The  basic  criterion  for  evaluating  gross  performance  trends 
from  the  single  sample  runs  was  to  determine  the  state  variable 
"error."  This  error  is  defined  as  the  difference  in  magnitude  between 
the  estimator  state  solution  and  the  solution  associated  with  the  non¬ 
noise  corrupted  "truth  model"  from  which  the  simulated  observations 
were  derived.  Results  which  show  an  error  significantly  greater  in 
magnitude  than  the  standard  deviation  of  the  estimator-computer  state 
covariance  values  were  used  for  case  selections  for  later  Monte 
Carlo  simulations.  For  the  sake  of  brevity  this  error  assessment  is 
presented  only  at  the  time  of  the  last  observation  along  the  reentry 
trajectory. 
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0.1.1.  Dynamics  Mismatch  Results 


A  number  of  truth  models  were  selected  to  examine  the  performance 
aspects  of  the  estimator  with  significant  mismatch  between  the  estima¬ 
tor  dynamics  model  and  the  functional  form  of  the  truth  model  dynamics. 
These  included  truth  models  with  dynamic  variations  of :  B  as  a  function 
of  Mach  number,  a  high  density  exponential  density  profile,  a  standard 
atmosphere  density  profile,  a  discrete  change  in  ballistic  parameter, 

3,  and  various  combinations  of  these  dynamics.  The  truth  model 
trajectories  and  the  dynamics  mismatch  profiles  of  the  variable  B  and 
atmospheric  density  profiles  are  contained  in  Appendix  B.l.  Appendix 
B.2.  contains  the  details  of  the  single  sample  estimator  runs  which 
are  summarized  below. 

The  initial  epoch  for  this  simulated,  decayed  satellite  tra¬ 
jectory  was  at  an  altitude  of  73.82  KM  with  a  0.5  degree  reentry  angle. 
Initial  positions  and  velocities  are  identical,  within  both  the 
estimator  and  the  truth  model.  Initial  estimator  conditions  for  x^  and 

xQ  vary  from  the  truth  model  only  if  that  particular  case  differs  from 
o 

the  estimator  dynamics. 

Most  single  sample  cases  included  exact  (non-noise  corrupted), 

truth  model  derived  observations  with  a  10~5  radian  one-sigma  weighting 

within  the  observation  covariance  matrix,  R,  , .  In  selected  cases, 

(n; 

corruptive  random  noise  was  added  to  the  simulated  observations  to 
examine  the  combined  effect  of  dynamics  mismatch  and  observation  noise 
on  the  solution  process.  With  the  observation  satellite  positioned  in 
a  synchronous  orbit  this  10-5  radian  noise  represents  a  position 
uncertainty  of  approximately  0.4  KM  at  the  points  along  the  reentry 
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trajectory.  The  noise  corrupted  simulations  were  further  examined  in 
the  subsequent  Monte  Carlo  analysis. 

The  single  sample  results  are  shown  in  Table  I  at  the  time  of 
the  last  observation,  T  =  330  seconds  trajectory  time,  illustrating: 
error,  estimator  covariance  standard  deviation  (one  sigma)  values, 
and  the  error/ (one  sigma)  ratio  for  each  pertinent  state  variable. 

This  single  sample,  error/(one  sigma)  ratio  will  be  referred  to  as  the 
"performance  ratio"  in  the  text.  The  one  sigma  value  is  from  the 
estimator-computed  state  covariance  matrix.  (Note  that  in  the  cases  of 
a  U.S.  1962  Standard  Density  truth  model  (ATM  62),  no  error  or  ratio 
comparison  was  made  for  the  exponential  scale  height,  x^,  since  the 
results  of  interest  are  the  position  and  velocity  values  determined  by 
using  an  exponential  atmospheric  density  model  while  processing  data 
derived  from  the  use  of  a  standard  atmosphere.) 

Under  Monte  Carlo  simulation  analysis,  significant  performance 
degradation  of  the  estimator  would  be  evidenced  by  a  growth  in  the 
bias/ (one  sigma)  ratio  much  beyond  one.  These  bias  values  represent 
the  mean  error  in  state  variables  over  the  n  Monte  Carlo  runs.  With 
the  qualifications  normally  extended  to  single  sample  assessments,  one 
can  draw  some  trend  information  by  review  of  these  results.  Concen¬ 
trating  on  the  performance  ratio,  results  of  the  various  truth  model 
cases  in  Table  I  show  the  following: 

1)  The  ratio  generally  increases  with  increasing  mismatch 
between  the  estimator  and  truth  model  dynamics,  as 
expected.  There  is  a  distinct  variation  in  which  states 
are  more  or  less  effected  which  is  related  to  the  particular 
truth  model  data  being  processed  by  the  estimator. 
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Single  Sample  Results:  Infinite  Memory,  Dynamic  Mismatch 
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2)  The  ratio  generally  increases  with  the  addition  of  random 


noise  to  the  observations.  This  is  not  a  true  indication 
of  the  estimator  performance  with  random  noise  on  the 
observations,  since  only  a  single  sample  trajectory  is 
used.  The  later  Monte  Carlo  results  present  a  more 
complete  illustration  of  the  estimator  performance  with 
noise  corrupted  observations. 

3)  The  ratio  generally  decreases  (improved  performance)  with 
one  observation  per  update  compared  to  two  observations. 

This  indicates  an  increasing  ability  of  the  deterministic 
estimator  dynamics  model  to  match  the  truth  model  over 
smaller  time  intervals  between  the  epoch  update  point  and 
one  set  of  observational  data.  This  provides  an  indication 
of  the  need  to  apply  some  form  of  dynamic  compensation  to 
the  deterministic  estimator  model.  Under  normal  circum¬ 
stances,  one  would  anticipate  an  improved  performance  with 
more  observations.  However,  the  deterministic  estimator 
dynamics  model  is  valid  over  a  very  limited  time  period. 

A)  The  large  ratio  and  error  values  of  the  step  change  in 
6  at  T  =  300  seconds  result  from  the  small  number  of 
observations  after  this  discrete  change  in  8  and  the  small 
state  covariance  values  resulting  from  having  already 
processed  a  large  number  of  observations.  The  dynamics 
model  was  then  excessively  weighted  in  the  update  process 
and  did  not  fully  accommodate  the  information  content  within 
the  final  few  observations.  Also,  the  small  number  of 
observations  occurring  after  this  step  change  were  not 
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sufficient  for  the  effects  of  the  transient  nature  of  this 
discrete  change  in  dynamics  to  settle  out  in  the  estimator 
solution.  Again,  this  provides  further  evidence  that  some 
dynamic  compensation  is  required  for  the  deterministic 
dynamics  model. 

Residual  monitoring,  which  is  often  used  to  help  determine  a 
pseudo-noise  matrix  or  fading  memory  parameter  selection,  was  considered 
to  address  the  limitations  in  the  estimator  dynamics  model.  As 
Morrison  (2)  and  Sorenson  &  Sacks  (3)  point  out:  observation  residual 
testing  is  most  valid  when:  i)  the  functional  form  of  the  model 
dynamics  is  sound,  and  ii)  sufficiently  large  numbers  of  residuals  are 
available  over  the  tine  span  of  local  dynamic  model  validity,  for  the 
statistical  analysis  of  the  residuals  to  be  valid.  Unfortunately,  as 
the  expanded  discussion  in  Chapter  III  will  show,  the  time  span  of 
local  model  validity  is  often  very  short.  The  implications  of  residual 
monitoring  with  a  fading  memory  estimator  are  discussed  in  more  detail 
in  Chapter  III. 
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D.1.2.  Observation  Noise  Variations 

With  the  same  initial  state  conditions  within  the  truth  model  and 

the  estimator,  a  sequence  of  simulated  runs  were  completed  to  analyze 

the  effect  of  variations  of  the  one  sigma  corruptive  noise  on  the 

observations.  The  same  dynamics  were  modeled  in  both  the  estimator 

and  the  truth  model.  Random  noise  with  the  same  initial  seed  was 

placed  upon  the  simulated  observations  using  three  different  levels  of 

observation  noise.  In  each  case  the  estimator  observation  covariance, 

R,  , ,  was  consistent  with  the  level  of  random  noise  on  the  observa- 
(n) 

tions.  The  resulting  error  and  performance  ratio  are  shown  in  Table  II 
at  T  =  160  seconds  trajectory  time,  the  time  of  the  final  observation. 

A  more  severe  nonlinear  trajectory  (i.e.,  lower  altitude,  higher 
atmospheric  density)  with  an  initial  altitude  of  53.73  KM  was  used  in 
this  comparison. 

The  trends  which  may  be  observed  from  a  review  of  the  Table  II 
results  include  the  following: 

1)  The  estimator  error  magnitude  improves  going  from  10“ 4  to 
10“ 5  radian  noise  levels. 

2)  Acceptable  estimator  solutions  in  terms  of  the  performance 
ratio  are  indicated  with  the  exact  dynamics  in  both  the  10-4 
and  10"5  radian  noise  cases.  The  performance  ratio  is  not 
significantly  greater  than  one.  The  onset  of  growth  in  this 
ratio  is  evident  with  a  10“5  radian  noise  level. 

3)  With  an  extremely  accurate  set  of  observations  (10  6  radian 
noise),  performance  degrades  when  compared  to  the  10“4  and 
10”5  radian  noise  cases.  This  is  likely  due  in  large  part 


to  the  overweighting  of  the  deterministic  estimator 


dynamics  model  late  in  the  trajectory  data  span,  due  to  a 
more  rapid  collapse  in  the  magnitudes  of  the  state 
covariance  terms  with  10”6  radian  observation  noise.  The 
linearity  check  of  Equation  44  was  consistently  satisfied 
throughout  the  trajectory  estimation.  This  provides  a 
reasonable  assurance  that  the  error  in  the  state  estimate 
was  not  primarily  due  to  higher  order  term  corruption. 
However,  model  compensation  or  a  fading  memory  may  be 
required  with  extremely  accurate  observations. 
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Single  Sample  Results:  Infinite  Memory,  Observer  Noise  Variations 
(1  obs  per  update  epoch  -  exact  dynamics) 


D.1.3.  Observer  Angle  Variations  Relative  to  the  Reentry  Trajectory 
With  the  same  dynamics  and  initial  state  conditions  within  the 
estimator  and  the  truth  model,  a  sequence  of  runs  was  made  to  compare 
the  effects  of  variations  in  the  orbital  location  of  the  observer 
relative  to  the  reentry  trajectory.  These  runs  were  obtained  with  an 
initial  trajectory  epoch  at  53.73  KM  and  may  be  compared  to  their 
respective  noise  level  counterparts  of  Table  II,  where  the  observer  was 
in  a  synchronous  orbit  with  an  initial  Right  Ascension  (RA)  of  +45°. 
Figure  3  shows  the  variations  in  the  initial  observer  orbital  positions 
relative  to  the  reentry  trajectory  case  considered. 
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FIGURE  3 


Table  III  shows  the  estimator  results  at  the  time  of  the  final 
rvation,  T  =  160  seconds.  A  review  of  these  results  shows: 

1)  There  are  no  drastic  variations  in  estimator  performance 
with  minor  variations  in  the  single  observer  location 
relative  to  the  reentry  trajectory. 

2)  A  subtle  relationship  between  individual  state  variables 

and  observer  location  indicates  some  potential  limitations 

of  observability  with  angle  only  observations.  This  is 

apparent  by  examining  the  performance  ratio  results  of  the 

10“ 5  noise  level  cases  at  initial  RA  of  22°  and  45°.  For 

example,  the  ratios  are  much  higher  for  the  x  and  y  position 

terms  and  the  Ep  term  with  the  22°  RA  case.  In  all  cases, 
o 

the  ratio  does  not  grow  beyond  one.  These  effects  are  more 

clearly  illustrated  in  Section  D.I.4.,  below.  As  the 

angular  separation  between  the  observer  subpoint  and  the 

reentry  satellite  becomes  smaller,  there  is  a  loss  in 

observability  available  from  the  angular  observations.  This 

is  not  an  unanticipated  result,  but  it  is  difficult  to 

examine  analytically  with  the  current  observer.  This  is  due 

to  the  coupling  of  the  geometric  partials  matrix  and  the 

state  transition  matrix  within  the  estimator  T,  .  term  (see 

(n) 

Equation  16).  The  articles  by  Sivazlian  and  Green  (21,22) 
lend  a  useful  insight  into  the  accuracy  considerations  of 
similar  angles  only  data  for  the  tracking  of  a  stationary 
target  from  multiple  observers.  Their  formulation  is  more 
easily  expressed  analytically.  They  show  a  general 
deterioration  in  the  accuracy  of  the  target  position 


estimates  as  a  function  in  the  angular  separation  of  the 


target  from  the  observers  and  the  angular  separation 
between  the  multiple  observers.  They  also  show  that  this 
function  relationship  loses  its  dependence  on  the  observer 
locations  if  the  observers  are  positioned  at  right  angles 
to  the  target.  When  extended  to  the  dynamic  case  of  the  cur¬ 
rent  observer,  similar  results  are  obtained.  Section  D.1.4. 
shows  the  loss  in  estimator  performance  as  the  reentry  data 
approaches  the  observer  subpoint.  Section  D.1.5.  shows  the 
improvement  in  estimator  performance  with  data  from  dual 
observers,  orthogonally  positioned  relative  to  the  reentry 
trajectory . 
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Single  Sample  Analysis:  Infinite  Memory,  Observer  Angle  Variations 
(one  observation  per  update  epoch) 


D.1.4.  Reentry  Trajectory  Location  Variations 


Because  of  the  trends  apparent  from  Table  III  with  observer 
angle  variations,  a  sequence  of  simulation  runs  was  completed  with 
variations  on  the  initial  trajectory  RA.  Observations  are  from  an 
observer  with  an  initial  RA  of  ±45°  in  each  case.  The  initial 
conditions  and  the  dynamics  models  were  matched  exactly  within  the 
estimator  and  the  truth  model.  A  10-5  random  noise  was  placed  upon 
each  set  of  observations.  A  lower  altitude  (49.86  KM),  high 
inclination  trajectory  (approximately  85°)  was  used  for  this 
comparison.  This  would  stress  the  estimator  performance  by 
i)  beginning  in  a  more  nonlinear  dynamics  region  of  higher  atmospheric 
density,  and  ii)  proceeding  on  such  a  high  inclination  trajectory 
would  result  in  smaller  changes  in  observation  azimuth  measurements, 
with  the  reentry  trajectory  case  more  closely  approaching  the 
observation  satellite  subpoint  or  nadir:  the  Earth  surface  location 
intersecting  with  the  vector  between  the  observation  satellite  and  the 
ECI  coordinate  system  origin  at  the  center  of  the  Earth.  These 
variations  in  reentry  trajectory  geometry  relative  to  the  observation 
satellite  are  shown  in  Figure  4. 

From  the  synchronous  altitude  of  the  observer,  the  figure  of 
the  Earth  appears  approximately  as  a  17°  solid  cone.  Therefore  in  the 
observer  coordinate  system,  the  angular  observation  of  any  reentry 
trajectory  can  at  most  vary  between  ±8.5°  relative  to  the  observation 
satellite  nadir.  The  results  of  this  "Earth  limb  to  satellite  nadir" 
variation  analysis  are  shown  in  Table  IV  at  a  trajectory  time  of  150 
seconds,  the  time  of  the  final  observation. 

A  review  of  the  Table  IV  results  reveals  the  following: 
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Single  Sample  Analysis:  Infinite  Memory,  Variations  in  Reentry  Trajectory  Position 
(one  observation  per  update  epoch,  exact  dynamics) 
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1)  There  is  a  general  deterioration  in  the  estimator  solutions 

in  terms  of  the  performance  ratio  as  the  trajectory  approaches 
very  close  to  the  observation  satellite  subpoint. 

2)  With  a  single  observation  satellite,  even  with  an  exact 
dynamics  model,  a  subtle  interrelationship  exists  between 
the  observer  geometry  and  the  trajectory  geometry  which  may 
affect  the  estimation  of  certain  states  more  than  others. 

This  may  be  seen  by  examining  the  differences  between  the 
RA  =  0°  and  45°  cases  for  the  individual  state  variables. 

As  the  observation  aspect  on  the  reentry  trajectory  varies, 
there  are  different  variations  in  both  error  and 
performance  ratio  for  the  individual  state  variables. 

A  rigorous  examination  of  these  variations  near  the  observer 

subpoint  is  an  extremely  complex  problem.  The  geometric  relationships 

are  coupled  to  the  particular  trajectory  via  the  product  = 

H(x  (t  ))  <Kt  ,  t  )  in  the  state  and  covariance  update  expressions 
mm  mm 

n  n 

(Equations  39  and  47).  A  review  of  Appendix  A.l.  and  Appendix  A. 2. 

shows  the  complexity  of  the  A(t)  matrix  used  to  generate  <f>  and  of  the 

geometry  matrix,  H(x  (t  )).  However,  one  can  obtain  some  insight  into 

mm 

n 

these  considerations  by  referring  to  the  observation  geometry  of 

Figure  1  and  the  R.A.  =  42°  results  of  Table  IV. 

Very  near  the  subpoint,  the  y'  and  observer  coordinates 

approach  zero.  As  Table  IV  shows,  the  estimator  solutions  deteriorate 

most  dramatically  in  the  three  state  position  terms  (x,  y  and  z)  and 

the  Bp  term.  Very  minor  errors  in  y'  and  z'  coordinates  couple  very 

directly  into  the  trajectory  position  coordinates  x,  y  and  z  through 

the  geometry  matrix,  H(x  (t  )).  Observability  of  the  velocity  states 

m  m 

n 
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is  obtained  principally  from  the  state  transition  matrix,  <j>,  within  the 

T^  matrix.  With  minimal  bias  in  the  velocity  and  the  density  scale 

height  components,  the  errors  in  position  generate  errors  in  the 

solution  for  gp  . 

o 

This  illustration  for  the  observer  subpoint  case  was  chosen  for 
examination  of  the  observability  limitations  from  a  single  observer 
providing  angular  data  in  this  region.  The  utilization  of  dual 
observers  may  be  advisable  to  improve  these  geometric  performance 
relationships  between  the  observation  satellite  and  reentry  trajectory. 
Since  in  the  general  satellite  decay  case,  no  control  exists  over  these 
geometry  relationships,  multiple  observation  satellites  may  be  required 
to  insure  acceptable  estimator  performance.  Examination  of  the 
contribution  of  a  dual  observation  satellite  set  of  measurements  will 
be  explored  in  the  next  set  of  single  sample  results. 


D.1.5.  Multiple  Observation  Satellites 

To  allow  direct  comparison  to  the  Section  D.1.3.  results,  the 
53.73  KM  initial  altitude  trajectory  was  used  to  complete  a  sequence  of 
single  sample  simulations.  These  runs  included  variations  in  the 
observer  initial  RA  and  observation  noise  levels,  with  data  from  two 
observers.  Table  V  contains  the  results  of  these  dual  observer  runs 
with  error  and  performance  ratio  values  at  T  =  160  seconds  trajectory 
time.  These  results  may  be  compared  to  the  single  observer  satellite 
results  of  Table  III.  Figure  5  shows  the  initial  orientation  of  the 
dual  observers  relative  to  the  reentry  trajectory.  This  observer 
geometry  depicts  a  situation  where  a  set  of  four  synchronous  observers 
are  positioned  in  orbit  to  provide  a  full  visibility  of  the  Earth. 

The  two  observers  shown  would  be  in  a  position  to  observe  the  selected 
reentry  trajectory. 

A  basic  consistency  with  the  previous  single  observer,  single 
sample  analysis  is  apparent,  with  some  additional  performance  benefits 
evident  using  data  from  two  observation  sources.  The  ±45°  initial  RA 
observers  provided  superior  estimator  performance  in  the  10” 4  radian 
observation  noise  cases.  The  10” 5  radian  noise  cases  for  ±45°  observer 
RA  indicate  a  subtle  combination  of  a  number  of  factors, 

1)  higher  observation  accuracy, 

2)  different  observation  geometry,  and 

3)  a  higher  weighting  of  the  estimator  dynamics  model  late 
in  the  observation  time  span. 

A  more  accelerated  collapse  in  the  magnitudes  of  the  estimator-computed 
state  covariance  entries  occurs  due  to  processing  twice  the  data  of  the 
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Single  Sample  Analysis:  Infinite  Memory,  Dual  Observer  Analysis 
t  dynamics,  2  observations  per  update  epoch  -  one  from  each  observer) 
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single  observer  case.  All  of  these  factors  begin  to  affect  the 
estimator  solutions  in  terms  of  the  performance  ratio  measure. 

Although  single  sample  runs  are  a  dangerous  basis  for  making  firm 
conclusions,  one  can  reasonably  expect  superior  estimator  performance 
from  the  dual  observer  solutions  —  unless  the  combination  of  higher 
data  volume  (more  observations)  and  more  accurate  data  combine  to 
generate  extremely  small  state  covariance  magnitudes  late  in  the 
trajectory.  In  this  event,  the  output  of  the  deterministic  estimator 
dynamics  model  will  be  overly  weighted,  thereby  ignoring  some  of  the 
information  content  of  the  later  observations.  Model  compensation  is 
then  required  to  correct  for  this  effect. 

The  combination  of  data  from  two  observation  sources  with  mis¬ 
matched  dynamics  and  a  fading  memory  estimator  formulation  will  be  more 
extensively  investigated  in  the  Monte  Carlo  numerical  results  of 
Chapter  III. 

D. 2.  Monte  Carlo  Analyses 

As  the  previous  single  sample  simulations  indicate,  the  major 
degradation  of  estimator  performance  occurs  in  the  presence  of  mismatch 
between  estimator  model  and  true  trajectory  dynamics.  While  other 
factors,  such  as  observation  geometry,  single  observation  satellite 
data,  and  highly  accurate  observations  ( 10~ 6  radians),  indicate  the 
performance  limits  of  the  estimator ;  the  dynamics  mismatch  cases 
clearly  contain  the  most  consistently  biased  solutions.  The  purpose  of 
the  Monte  Carlo  runs  presented  in  this  section  Is  to  demonstrate  a 
statistically  sound  quantification  of  the  performance  impacts  with  an 
exact  match  and  a  mismatch  between  the  estimator  and  the  truth  model 


dynamics . 


A  baseline  for  the  Monte  Carlo  analysis  was  established  as 


follows : 

1)  A  high  altitude  (73.82  KM)  initial  epoch  trajectory  was  used 
with  state  covariance  matrix  values  consistent  with  position 
and  velocity  uncertainties  propagated  from  a  final  predecay, 
orbital  solution.  Covariance  values  for  Bp  (x.,)  and  Q  (xQ) 
were  consistent  with  either  nominal  uncertainties  for  the 
matched  dynamics  cases  or  with  the  initial  conditions 
represented  by  the  dynamics  mismatch  between  the  estimator 
and  the  truth  model. 

2)  The  values  of  the  initial  state  variables  were  randomly  varied 
in  accordance  with  the  initial  state  covariance. 

3)  The  simulated  data  consisted  of  a  set  of  33  angular  observa¬ 
tions  for  each  case,  at  10  second  intervals,  randomly 
corrupted  with  either  10-1*  or  10-5  radian  observation  noise. 
All  data  came  from  a  single,  synchronous  observation 
satellite  with  an  initial  RA  =  +45°. 

4)  Each  series  of  runs  retained  identical  seeds  on  the  random 
initial  state  variables  variations  and  the  observational  data 
for  valid  comparison  between  cases. 

5)  Cases  involving  exact  knowledge  of  the  dynamics  were  completed 
with  two  observations  sets  per  update  epoch.  The  dynamics 
mismatch  case  was  completed  with  one  observation  set  per 
update  epoch  to  insure  maximum  local  validity  of  the 
estimator  during  state  update  and  for  propagation  between 
observations.  Additional  simulations  showed  that  up  to  5  or 
10  observations  could  be  processed  for  the  initial  epoch 
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update  in  the  near  circular  orbit  conditions  early  in  the 


trajectory.  As  the  maximum  aerodynamic  deceleration  regions 
are  reached,  this  time  span  over  which  the  estimator  update  is 
valid  is  very  limited.  For  ease  of  implementation  on  the 
digital  computer,  this  update  time  span  was  restricted  for  the 
worst  case  along  the  entire  trajectory.  An  attempt  was  made 
to  dynamically  limit  the  number  of  observations  used  to  update 
the  trajectory  epoch  points  by  an  application  of  the  lineariza¬ 
tion  validity  check  of  Equation  44.  This  proved  difficult  due 
to  variations  which  occurred  in  the  numbers  of  admissible 
observations  changing  as  a  function  of  the  iteration  number  in 
a  given  update.  As  the  differential  corrector  converged, 
different  numbers  of  observations  could  be  used  without 
violating  the  criteria  of  the  linearization  check. 

6)  The  dynamics  mismatch  results  included  simulated  observations 
from  a  truth  model  with  8  =  f(Mach  no.)  and  atmospheric 
density  from  the  U.S.  1962  Standard  Atmosphere.  These 
represented  a  first  step  towards  processing  the  functional 
variations  representative  of  "true  reentry  dynamics"  relative 
to  the  more  simplified  estimator  model. 

7)  The  sequence  of  30  samples  was  selected.  This  number  was 
chosen  because  a  negligible  change  occurred  in  the  mean  bias 
values  of  the  state  estimates  with  higher  numbers  of 
replications.  An  example  will  illustrate  the  relatively 
insignificant  change  in  the  bias  magnitude  while  processing 
the  Monte  Carlo  runs.  The  next  page  shows  the  infinite 
memory  bias  magnitudes  for  the  dynamics  mismatch  Monte  Carlo 
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results  on  runs  number  27  through  30.  These  may  be  compared 


to  the  mean  value  and  standard  deviations  computed  in  the 
Monte  Carlo  analysis  to  see  there  is  a  rather  minor  change  in 
mean  bias  value  with  a  further  increase  in  numbers  of 
replications:  minor  as  compared  to  the  mean  and  standard 
deviation  of  the  Monte  Carlo  results. 


10"5  Mismatch  Dynamics  CASE  3  INF 
Monte  Carlo  Reduction 


Bias  after 

Bias  after 

Bias  after 

#27 

#28 

#29 

X 

3.2410 

3.1238 

3.08196 

X 

1.810  X  10'5 

1.7823 

X 

10- 5 

1.8093  X  10“ 5 

y 

3.3603 

3.2378 

3.1894 

y 

2.274  X  10"5 

2.2108 

X 

10“ 5 

2.3513  X  10'5 

z 

.2514 

.2446 

.2457 

z 

2.6597  X  10"6 

2.5606 

X 

10“6 

2.5152  X  lO"6 

ep0 

2.5596  X  10- 3 

2.5095 

X 

10-3 

2.4766  X  10"3 

Q 

4.0485  X  lO"3 

4.0471  X 

10“ 3 

4.1245  X  10'3 

Bias  after 
#30 

Mean 

Estimator 
Standard  Deviation 

X 

3.0037 

994.67 

2.1958 

♦ 

X 

1.7514  X  10“5 

-3.6205 

.0068142 

y 

3.1155 

2053.8 

2.1960 

• 

y 

2.4319  X  lO'5 

4.2434 

.004735 

z 

.2375 

5961.8 

.6064 

z 

2.42996  X  10'6 

-1.4195 

.0017136 

epo 

2.4688  X  10~3 

.50675 

.057406 

Q 

4.0343  X  10“3 

6.9917 

.063271 

Trajectory  time  T  =  330  seconds 
at  last  observation 
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Specific  details  on  the  estimator  model  initial  conditions  are 


shown  in  Table  VI  below. 


Table  VI 

Estimator  Initial  Epoch  Conditions 
(Reference  Trajectory) 


Altitude:  78.82  KM 
RA:  0° 

Declination:  68.1° 

Inclination:  10.9° 

x  -  2396  KM  8p  -  .49* 

o 

X  -  -3.905  KM/sec  Q  -  7.0031  KM 

y  -  0  KM 
y  -  6.67  KM/sec 
z  -  5965  KM 
z  -  1.49  KM/sec 


★ 

Equivalent  to  6 


4  X  10“10  KM2/KG 


Complete  details  on  the  6  =  f (Mach  no.)  and  U.S.  1962  Standard 
Atmosphere  are  included  in  Appendix  B.l. 
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The  specification  of  conditions  for  convergence  of  each  iteration 


of  the  estimator  at  the  epoch  updates  was  consistent  with  existing 

orbital  applications  (5,13).  Values  were  chosen  to  preclude  generating 

excessive  numbers  of  estimator  iterations  with  insignificant  changes  in 

the  6x(t  ).  As  mentioned  earlier  in  Section  A.,  the  cj  and  zz  were 
m 

empirically  determined  such  that  the  change  in  individual  state  solu¬ 
tions  was  small  (less  than  approximately  .01  times  the  standard  devia¬ 
tion)  on  the  final,  converged  iteration.  The  estimator  convergence 
criteria  were  specified  as  follows: 

-  Relative  converge  criterion,  (ei  Equation  34)  equal  to: 


.20 

Matched  dynamics 

10"  4 

observation 

noise 

(CASE 

11NF) 

.15 

Matched  dynamics 

10"5 

observation 

noise 

(CASE 

2  INF) 

.15 

Mismatched 

dynamics 

10"5 

observation 

noise 

(CASE 

3  INF) 

-  Absolute  convergence  criterion  of  (e2  Equation  35) : 


.15 

X 

io-1* 

Matched 

dynamics 

10'4 

observation 

noise 

.10 

X 

10'5 

Matched 

dynamics 

10"5 

observation 

noise 

.10 

X 

10"  5 

Mismatched  dynamics 

10'5 

observation 

noise 

The  total  change  in  the  state  variables  5x(t  ).  was  then  examined 

mi 

relative  to  their  a  priori  standard  deviations,  /  S  .  .  This 

m,m-l .  . 

li 

criteria  was  violated  on  a  significant  number  of  updates  only  for  the 
mismatched  dynamics  case. 

Each  of  the  three  Monte  Carlo  cases  are  denoted  by  case  number 
with  their  respective  initial  state  covariance  values  in  Table  VII, 
below.  This  table  shows  a  number  for  each  infinite  memory  case  (e.g., 
1INF),  the  observation  noise  standard  deviation,  the  number  of  observa¬ 


tions  per  epoch  update,  and  the  standard  deviations  from  the  estimator 


Table  VII 


Monte  Carlo  Baselines 


Initial  one  sigma 

do) 

covariance 

Case  No. 

Assumptions 

value 

1INF 

10-1*  obs  noise 

X 

4.511  KM 

2  obs/epoch  update 
matched  dynamics 

X 

.1  KM/sec 

y 

9.268  KM 

y 

.1  KM/sec 

Total 

position 

la: 

11.28  KM  (mostly  intrack) 

z 

4.575  KM 

.173  KM/sec 

z 

.1  KM/sec 

Total 

velocity 

la: 

.1  (20%  from 

O 

mean) 

Q 

.1  (16%  from 

mean 

density  p) 

2INF 

10“5  obs  noise 

X 

1.427  KM 

2  obs/epoch  update 
matched  dynamics 

• 

X 

.0316  KM/sec 

y 

2.931  KM 

• 

y 

.0316  KM/sec 

Total 

position 

la : 

3.566  KM  (mostly  intrack) 

z 

1.447  KM 

.0548  KM/sec 

• 

z 

.0316  KM/sec 

Total 

velocity 

la: 

6po 

.1 

Q 

.1 

3INF 

10" 5  obs  noise 

X 

1.427  KM 

1  obs/epoch  update 
mismatched  dynamics 

X 

.0316  KM/sec 

y 

2.931  KM 

y 

.0316  KM/sec 

Total 

position 

la: 

3.566  KM  (mostly  intrack) 

z 

1.447  KM 

.0548  KM/sec 

• 

z 

.0316  KM/sec 

Total 

velocity 

la: 

ep 

.316  (65%  from 

o 

mean) 

Q 


.1 


state  covariance  matrix  for  the  initial  state  conditions.  The  values 
that  these  standard  deviations  represent  for  3pq  and  Q  are  shown  as  a 
percent  deviation  from  the  mean  values  for  and  atmospheric  density, 

P. 


The  results  of  the  Monte  Carlo  analysis  are  shown  in  Figures 
6-17.  Total  position  and  velocity  data  are  shown  in  the  Monte  Carlo 
figures  for  ease  of  viewing  the  eight  state  system.  The  definitions 
for  the  presentation  and  discussion  of  the  Monte  Carlo  results  are  shown 
below  using  the  position  as  an  example: 


Total  position  bias: 

/ 

n 


BIAS  = 


n 


2  f(xt  -  X.)'  +  (yt  -  y.r  +  (zt  +  z^2] 


2-1^1 


(74) 


where:  x^  =  true  x  position  term 

=  single  estimator  run  estimate  for  x  position  term 
n  =30  replications 
1  sigma  (standard  deviation)  measures: 

1)  From  the  estimator  computed  covariance: 

1  sigma  est.  2  [o^2  +  oy2  +  o^2]2 

These  values  are  derived  from  the  average  of  the  state 

covariance  magnitudes  over  the  30  estimator  runs.  They  come 

from  the  updated  state  covariance,  S  ,  to  maintain 

m,m 

consistency  with  the  updated  position  and  velocity  values 
used  to  present  these  results  (Figures  6-17). 

2)  RSS  =  the  root  sum  square  of  the  Monte  Carlo  derived  mean 


square  errors  from  the  true  value. 


(75) 


RSS  =  [  —V  E  (x  -xj't  —V  Z  (y  -  y.)' 
n-1  .  ,  t  i '  n-1  .  ,  Jt  Ji 

i=l  1=1 

,  n  , 

+  ^r  =  -  zi)2i" 

1=1 

3)  1  sigma  about  mean  solution  =  RSS  of  the  variance  about  the 
mean  solution. 


=  [  — Y  I  (x  -x,)2^-^  £  (y  -  y.)‘ 

n-1  .  .  m  i  n-1  .  ,  m  l 
i=l  i=l 


(76) 


,  n  -  i 

+  I  (z  -  z.)2]"1 

n-1  .  '  m  i  1 

i=l 


where:  x  =  mean  estimator  solution  for  x  position  term  for 
m 

n  runs 

Similar  expressions  were  developed  for  the  total  velocity  magnitude 
terms . 

Two  principal  measures  of  merit  were  used  to  assess  the  Monte 
Carlo  results.  In  the  figures  labeled  "ESTIMATOR  PERFORMANCE",  the 
ratio  of  RSS  of  the  mean  square  error  about  the  true  solution  to  ONE 
SIGMA  (average  estimator-computed  standard  deviation  value)  are  shown. 
This  ratio  gives  an  indication  of  the  validity  of  the  estimator 
covariance  matrix  values.  For  acceptable  estimator  performance  the 
ratio  should  remain  close  to  one. 

In  the  set  of  figures  labelled  "MONTE  CARLO  RESULTS",  the  basic 
objective  was  to  illustrate  the  mean  bias  magnitude  (over  30  samples) 
relative  to  the  estimator-computed  standard  deviation  values.  Two 
additional  measures  of  the  estimator  performance  are  also  displayed  to 
assist  in  assessment  of: 
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1)  The  validity  of  the  estimator-computer  state  covariance,  and 

2)  The  "apparent"  divergence  of  the  estimator  solution  in  the 
presence  of  mismatched  dynamics  models.  This  apparent 
divergence  is  defined  as  a  bounded  divergence  in  the  state 
solution.  This  was  the  -'haracter  of  the  growth  in  the 
position  and  velocity  biases  in  the  mismatched  dynamics  case 
3INF.  A  "true"  divergence,  by  contrast,  would  illustrate  an 
unbounded  growth  in  the  solution  bias.  By  strict  definition, 
a  true  divergence  would  require  some  specific  time  to  reach 
an  unbounded  condition.  Practically  speaking  for  this  appli¬ 
cation,  a  true  divergence  would  generate  a  skip  trajectory 
which  does  not  reenter  the  Earth  atmosphere. 

The  additional  measures  included  the  second  order  statistics 
(standard  deviation  values)  derived  from  the  30  Monte  Carlo  replicates. 
One  is  from  the  mean  square  errors  about  the  true  solution  (Equation 
75)  and  the  other  is  from  the  variance  about  the  30  sample  mean  solution 
(from  Equation  76).  With  the  growth  in  bias  in  the  estimator  solution, 
these  two  measures  of  the  estimator  performance  will  diverge.  The 
standard  deviation  about  the  true  solution  (Equation  75)  should 
increase  with  the  bias  in  the  estimator  solution.  If  the  standard 
deviation  about  the  mean  solution  (Equation  76)  maintains  levels 
consistent  with  the  estimator-computed  values  (Equation  51),  the 
estimator-computed  variance  will  be  indicative  of  the  real  random  error 
in  the  estimator  solution.  One  then  must  address  the  impact  of  the 
systematic  error  resulting  from  the  bias  in  the  estimator-computed 
state  variable  solutions.  Classical  methods  include  tuning  the 


estimator  so  that  its  computed  covariance  matches  actual  mean  square 


value  errors,  incorporating  "bias"  correction  terms  from  higher  order 
filters,  employing  full  scale  higher  order  filters,  etc.  (23,24,25) 


Examination  of  Figures  6-9  shows  acceptable  estimator  statistics 
are  available  from  the  estimator  covariance  matrix  with  matched 
estimator  and  truth  model  dynamics.  The  vertical  scale  on  these  figures 
is  exaggerated  to  a  ratio  of  10  for  comparison  to  the  mismatch  dynamics 
results  (Figures  14-15)  shown  later.  Similarly,  Figures  10-13  show  the 
position  and  velocity  bias  magnitudes  well  below  the  estimator  one 
sigma  levels.  The  results  also  substantiate  that  the  estimator  variance 
data  follows  a  trend  consistent  with  the  variance  data  derived  from  the 
Monte  Carlo  results  (Equations  75  and  76) . 

The  position  variance  results  show  a  growth  in  magnitude  as  the 
trajectory  approaches  the  170  second  point.  This  is  near  the  region  of 
the  decayed  trajectory  departing  the  "near  circular”  orbit  conditions 
and  into  a  committed  reentry.  As  the  maximum  deceleration  region  is 
passed  and  the  denser  lower  altitudes  are  reached,  the  position 
uncertainties  decrease  in  magnitude  (Figures  10,12). 

The  velocity  variance  values  show  an  almost  monotonic  decrease 
with  increasing  atmospheric  density  (Figures  11,13).  The  variations  on 
the  Monte  Carlo  derived  second  order  statistics  used  to  validate  the 
estimator  covariance  data  are  due  to  the  30  sample  replication  set  size. 
While  bias  magnitude  was  extremely  well  represented  by  30  samples,  at 
some  specific  time  points  the  Monte  Carlo  derived  variances  do  not 
follow  the  smooth  trend  of  the  estimator-computed  variance  data.  One 
factor  which  may  contribute  to  the  smoother  variations  in  estimator 
variance  data  along  the  trajectory  is  the  fact  that  the  values  shown 
are  average  standard  deviation  levels  over  the  30  runs,  rather  than 
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standard  deviations  computed  from  a  single  run.  In  any  event,  the 
trends  are  consistent  among  all  three  methods. 


ESTIMATOR  PERFORMANCE 


FIGURE  6 


FIGURE  7 


H 


UN9IS  3J0/SSH 


6  3dfK)U 


riGURE  10 


FIGURE  11 


FIGURE  12 


MONTE  CARI.O  RESULTS  -  CASE 


Clearly,  the  results  shown  in  Figures  6-13  demonstrate  valid  and 
acceptable  estimator  performance  with  matched  dynamics  at  two  different 
levels  of  observation  noise  and  initial  state  covariance  magnitudes. 
Unfortunately,  real  decay  trajectories  are  not  so  precisely  modeled. 

The  results  of  the  third  set  of  Monte  Carlo  runs  using  the  dynamic  mis¬ 
match  truth  model  are  shown  in  Figures  14-17. 

As  expected,  an  apparent  divergence  of  the  estimator  solution  is 
evident  in  both  the  "ESTIMATOR  PERFORMANCE"  ratio  of  RSS/(ONE  SIGMA) 
results  and  the  "MONTE  CARLO  RESULTS".  Both  bias  and  true  solution 
(Equation  75)  derived  standard  deviation  magnitudes  depart  from  the 
estimator  computed  variance  levels.  The  estimator  standard  deviations 
closely  follow  the  Monte  Carlo  derived  second  order  statistics  (Equation 
76)  about  the  30  sample  mean  solution.  These  results  show  a  clearly 
biased  solution  relative  to  the  estimator  indicated  statistics.  The 
primary  mechanism  driving  these  results  is  the  unmatched  dynamics 
between  the  estimator  and  truth  model.  The  solution  begins  to  diverge 
as  the  trajectory  bends  from  the  "near  circular  orbit"  conditions. 

Early  in  the  trajectory,  the  variations  in  6  were  relatively  small  and 
the  exponential  density  model  of  the  estimator  could  locally  represent 
the  U.S.  1962  Standard  Atmosphere  density  profile,  since  most  occur 
within  a  single,  non-isothermal  layer  of  the  standard  atmosphere. 

As  the  reentry  proceeded,  the  exponential  atmosphere  model  could 
not  completely  locally  accommodate  the  changing  density  gradients  of  the 
standard  atmosphere  which  contains  altitude  layers  of  positive  and 
negative  thermal  lapse  rate,  with  an  intervening  isothermal  altitude 
band.  Also,  the  more  significant  variations  in  the  truth  model 
B  =  f(Mach  no.)  occur  later  in  the  data  span.  Deceleration  in  the 
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trajectory  causes  the  velocity  to  decrease  from  high  hypersonic  values 
to  lower  supersonic  values  (Appendix  B.I.). 

An  additional  factor  which  substantially  influences  this  bias 
effect  is  the  inability  of  an  infinite  memory  estimator  with 
deterministic  dynamics  to  adjust  to  large  changes  late  in  the  data  span. 
This  is  a  natural  consequence  of  both  the  deterministic  dynamics  and  of 
reductions  in  the  estimator  covariance  terms  (particularly  velocity) 
with  the  accumulation  of  more  and  more  observations.  The  later  observa¬ 
tions  become  increasingly  less  influential  on  the  state  estimate,  with 
increased  weighting  on  the  output  of  the  estimator  dynamics  model. 
Unfortunately,  this  phenomenon  occurs  in  a  region  of  reentry  where  the 
most  marked  difference  between  the  truth  model  and  the  estimator 
dynamics  exists. 

Ideally,  one  would  desire  to  implement  a  pseudo-noise  compensa¬ 
tion  to  the  deterministic  model  dynamics.  The  magnitudes  of  the  noise 
would  then  be  tuned  to  minimize  the  bias  in  the  state  solution,  by 
matching  the  variances  in  the  estimator-computed  covariance  to  the  true 
mean  square  error.  Unfortunately,  this  procedure  has  not  met  with  much 
success  without  significant  amounts  of  empirical  trajectory  observations 
to  assist  in  the  tuning  process  for  the  noise  strength  (5).  The 
discussions  of  Chapter  III  show  the  model  compensation  methods  considered 
to  address  this  problem. 


ESTIMATOR  PERFORMANCE 


FIGURE  14 


ESTIMATOR  PERFORMANCE 
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FIGURE  10 


MONTE  CARLO  RESULTS  -  CASE 


D . 3 .  Special  Numerical  Investigations 

Several  special  simulation  runs  were  completed  to  illustrate 
important  influences  on  the  performance  aspects  of  the  basic  estimator 
formulation  in  addition  to  the  dynamics  mismatch  of  Figures  14-17.  One 
comparison  shows  the  difference  between  the  eight  dimensional  formula¬ 
tion  with  an  exponential  atmosphere  and  a  seven  dimensional  system, 
where  atmospheric  density  is  obtained  from  a  1962  U.S.  Standard 
Atmosphere.  Another  investigation  centers  about  the  proposed  eight 
variable  system.  This  illustrates  the  observability  considerations  and 
ill-conditioned  information  matrix  of  the  estimator  with  the  use  of 
angle  only  observations  from  a  single  observation  satellite.  Lastly, 
one  must  explore  the  validity  of  propagating  the  state  covariance 
matrix  from  the  last  epoch  to  impact  via  linear  techniques  (Equation 
46). 

An  illustration  of  how  the  more  complex  standard  atmosphere 
equations  contribute  to  corrupting  the  estimator  solutions  is  shown  in 
Table  VIII  showing  the  results  of  a  single  trajectory  estimation.  Each 
column  reflects  the  error  between  integration  of  the  nonlinear  state 
dynamics  equation  and  a  linear  propagation  of  the  5x(tm)  epoch 
correction  by  use  of  the  state  transition  matrix.  These  contrast  the 
results  of  the  eight  and  seven  dimensional  formulations  in  propagation 
of  a  typical  Sx  (x  position  term)  from  an  epoch  at  73.45  KM  altitude. 
Both  of  these  infinite  memory  estimators  had  matching  dynamics  between 
the  estimator  and  truth  model.  They  each  were  used  to  process  10 
observations  in  a  single  epoch  update.  The  resulting  state  correction 
was  propagated  forward  by: 
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1)  Propagate  forward  from  the  epoch  time,  t  ,  both  the  a  priori 

m 

and  updated  state  estimates,  x  (t  ,)  and  x  ,.(t  ,  ,),  by 

m  m+I  m+i  m+i 

integration  of  the  state  dynamics  equation.  One  can  find  the 
resulting  <5x(t)  by  differencing  these  two  solutions. 

2)  Propagate  the  <Sx(tm)  via  linear  methods,  i.e.,  6x(t)  = 

4>(t,t  )  <Sx(t  ).  This  is  the  form  of  the  incorporation 
m  m 

within  Che  state  update.  Equation  29. 

The  difference  in  these  two  propagations  gives  a  measure  of  the  error  in 
the  linear  approximation  of  the  estimator. 

Table  VIII 

ERROR  in  PROPAGATION 
(x  position  term) 


Trajectory  Time 

Observation 

ERROR 

8  dimensional  7  dimensional 

(sec)  from  Epoch 

Number 

system 

system 

10 

1 

.00002 

.0410 

20 

2 

.00006 

.0220 

30 

3 

.00013 

.0150 

40 

4 

.00022 

.0110 

50 

5 

.00033 

.0091 

60 

6 

.00046 

.0074 

70 

7 

.00062 

.0071 

80 

8 

.00080 

.0052 

90 

9 

.00099 

.0043 

100 

10 

.00122 

.0036 

The  errors  in 

the  eight  dimensional  formulation 

increase 

monotonically  with  time,  illustrating 

the  accumulation 

of  error  over  the 

100  seconds  from  the 

update  epoch.  This  underscores  the  previous 

recommendations  of  using  a  linearity  check  between  the  state  and 
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observations  changes  during  the  update  (Equation  44)  to  limit  the 
update  data  span.  By  limiting  this  time  span,  a  small  number  of 
observations  are  processed  at  each  epoch  to  minimize  the  effects  of 
higher  order  term  corruption  on  the  state  estimate. 

The  error  in  the  seven  dimensional  formulation  starts  at  a 
larger  magnitude  and  decreases  with  time.  This  is  a  reflection  of  the 
more  complex,  nonlinear  standard  atmosphere  density  equations  (49  and 
50)  in  the  various  altitude  layers  covered  by  this  propagation.  Recall 
that  the  eight  dimensional  algorithm  estimates  Q  as  in  Equation  51,  but 
the  seven  dimensional  estimator  does  not.  At  the  10  second  point,  the 
thermal  lapse  rate  is  a  -4.0,  near  the  center  of  the  61-79  KM  layer  of 
the  standard  atmosphere  (Appendix  B.I.).  The  remaining  propagation 
cuts  across  the  61-52  KM  layer  (-2.0  lapse  rate),  and  into  the 
isothermal  layer  of  47-53  KM  altitude.  The  errors  are  larger  in  the 
higher  thermal  lapse  rate  regions  and  are  higher  throughout  when  com¬ 
pared  to  the  eight  dimensional  formulation.  Since  the  magnitude  of  the 
error  can  be  large  in  a  given  non-isothermal  atmosphere  layer  the  use 
of  the  standard  atmosphere  approach  would  yield  significant  error  in 
the  state  solution,  even  with  a  limit  on  the  update  data  span.  There¬ 
fore,  one  must  consider  alternate  estimator  formulations  (e.g., 
statistical  linearization  approaches)  or  a  much  higher  observation 
data  rate  to  use  this  standard  density  model  formulation. 

The  system  observability  limitations  and  information  matrix  ill- 
conditioning  are  highlighted  by  examining  the  eight  dimensional 
formulation  using  different  numbers  of  observations  for  an  epoch  update. 
In  the  absence  of  a  priori  covariance  matrix  information,  a  minimum 
number  of  observations  must  be  processed  at  an  epoch  if  the  information 
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T 

matrix,  (T,  .  R.  S-1T,  .  )»  is  to  be  inverted  successfully,  and  thereby 
in;  in;  in; 

represent  a  valid  covariance  matrix.  Table  IX  shows  a  typical  example 

of  an  eigenvalue  analysis  on  the  information  matrix,  contrasting  a 

single  observation  update  with  a  10  observation  update  at  the  same 

epoch,  t  . 

m 

Table  IX 


Information  Matrix  Eigenvalues 


Single  Observation 
Update 


6.5225  X  10"15 
3.2239  X  10"1V 
7.4055  X  10"12 
7.9587  X  10-11 
6.1812  X  10"3 
1.2684  X  10“2 


10  Observation 
Update 

3.6973  X  10"6 
7.4146  X  10~6 
.2186 
12.6819 
176.5283 
2247.3215 
231252.1497 
10098973.1735 


(  Values  are  approximately  zero  due  to  computer 
roundoff  in  this  application.) 


The  single  observation  update  is  numerically  non-positive 
definite.  By  contrast  the  10  observation  update  is  numerically 
positive  definite.  In  either  case  the  information  matrix  is  highly 
ill-conditioned.  This  is  a  product  of  the  angle  only  observations 
from  a  single  synchronous  observation  satellite.  The  last  two 
eigenvalues  of  the  single  observation  update  reflect  the  observability 
of  the  system  available  from  a  single  two  dimensional,  angular 
measurement.  The  remaining  six  states  are  practically  not  observable. 
Even  the  10  observation  update  shows  a  13  order  of  magnitude  difference 
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in  the  size  of  the  eigenvalues.  The  two  smallest  eigenvalues  (of  order 
10“6)  result  from  the  inability  of  the  angular  measurements  to  observe 
the  position  and  velocity  along  a  line  of  sight  from  the  observer  to 
the  reentry  satellite. 

The  results  of  Table  LX  also  illustrate  that  a  certain  minimum 
number  of  observations  must  be  included  in  the  state  update  at  the 
initial  epoch  in  the  absence  of  a  priori  covariance  data.  Fortunately, 
the  early  observations  occur  over  a  region  with  a  small  reentry  angle 
(near  orbital  conditions)  where  processing  a  number  of  observations  in 
a  single  update  does  not  violate  the  linearization  validity  check. 

This  number  of  observations  used  for  the  initial  epoch  update  will 
depend  upon  the  particular  reentry  case  and  observation  geometry 
being  used.  However,  the  poor  observability  from  the  single  observer, 
angles  only  data  almost  necessitates  the  existence  of  a  priori  covari¬ 
ance  data. 

The  final  special  numerical  investigation  of  the  basic  estimator 
concerns  the  comparison  of  the  linear  propagation  of  the  last  epoch 
covariance  to  impact  (Equation  46)  or  a  Monte  Carlo  derived  impact 
covariance  matrix.  Use  of  a  sequential  estimator,  processing  small 
numbers  of  observations  per  epoch  update,  was  necessary  to  avoid 
accumulated  nonlinear  effects  from  corrupting  the  state  update.  The 
propagation  of  the  state  covariance  matrix  often  occurs  over  larger 
time  spaces  between  the  final  epoch  and  impact.  While  the  state 
update  can  be  influenced  by  the  joint  nonlinearities  within  the 
observation  geometry  and  the  model  dynamics,  the  covariance  propagation 
is  subject  only  to  dynamic  nonlinear  effects  on  the  state  transition 
matrix,  <J>(t,t  )•  If  the  nonlinearities  within  this  nonobservable 
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portion  of  the  trajectory  between  the  last  observation  update  epoch  and 
Earth  impact  significantly  affect  the  validity  of  j>(t,tm),  Equation  46 
cannot  be  used  to  propagate  the  state  covariance  over  this  region.  The 
alternative  of  a  Monte  Carlo  derived  impact  covariance  is  then 
available . 

The  Monte  Carlo  covariance  is  derived  from  a  90  replicate 
integration  of  the  state  dynamics  Equation  45  with  the  second  order 
statistics  accumulated  by  use  of  Equation  76.  A  random  sample  of  90 
variations  on  the  initial  conditions  from  the  final  epoch  covariance 
was  chosen  to  provide  a  more  true  representation  of  the  second  order 
statistics  over  the  large  time  spans  of  the  final  non-observable 
portion  of  the  trajectory.  Table  X  shows  the  comparison  of  these  two 
propagation  options,  by  presenting  the  standard  deviations  in  the 
state  variables.  These  were  obtained  by  propagating  from  the  last 
epoch  time  (T=324  seconds)  to  impact,  at  T=387.77  seconds.  This 
comparison  used  the  results  of  a  typical  estimator  run  with  10-5 
radian  one  sigma  observation  noise  and  exact  dynamics. 

The  results  of  Table  X  support  the  use  of  a  Monte  Carlo 
propagation  to  impact.  With  the  linear  propagation,  the  standard 
deviations  depart  from  the  Monte  Carlo  derived  values.  This  occurs 
initially  in  the  velocity  terms,  then  also  in  the  position  terms  as 
the  atmospheric  density  increases  and  the  velocity  magnitudes 
decrease.  This  Monte  Carlo  propagation  typically  requires  only  1.5 
to  1.7  times  the  computer  execution  time  over  the  linear  propagation. 

It  preserves  the  integrity  of  the  estimator  statistics  over  the  final 
nonlinear  region  of  the  trajectory  propagation. 


State  Covariance  Matrix  Propagation 


D . 4 .  Summary 

To  summarize  the  results  of  the  basic  estimator  findings,  key 
points  of  the  various  performance  investigations  are  detailed  in 
Table  XI. 


Table  XI 

Infinite  Memory  Estimator  Performance 

1.  Increasing  the  mismatch  between  the  estimator  and  the  truth  model 
dynamics  yields  the  most  dramatic  "bias"  in  the  estimator  state 
solutions . 

2.  With  deterministic  dynamics,  observation  noise  levels  of  less 
than  10-5  radian  (one  sigma)  begin  to  induce  significant  errors 
into  the  estimator  solution,  relative  to  the  standard  deviations 
of  the  estimator-computed  state  covariance  matrix. 

3.  With  identical  dynamics  in  the  estimator  and  truth  model: 

a.  Estimator  performance  is  acceptable.  The  bias  in  the  solution 
is  well  below  the  magnitudes  of  the  standard  deviations  of  the 
estimator-computed  state  covariance  matrix. 

b.  By  limiting  the  number  of  observations  used  for  an  update,  the 
effect  of  higher  order  term  corruption  on  the  estimator 
solutions  can  be  minimized. 

c.  Estimator  performance  deteriorates  as  the  reentry  trajectory 
approaches  the  subpoint  of  a  single  observation  satellite. 

d.  The  eight  dimensional  (exponential  atmosphere)  formulation  is 
superior  to  the  seven  dimensional  (standard  atmosphere) 
formulation,  in  the  current  estimator  using  the  Taylor's 
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series  linearization  approach  with  a  10  second  observation 


interval.  A  statistical  linearization  approach  may  be  more 
appropriate  for  the  seven  dimensional  system. 

e.  A  marginal  observability  with  the  angle  only  observations  from 
a  single  source  requires  existence  of  a  priori  covariance 
information  prior  to  an  epoch  update.  In  the  very  early, 
shallow  reentry  angle  portion  of  reentry,  a  state  covariance 
matrix  may  be  developed  from  batch  processing  a  number  of 
observations  at  the  initial  epoch.  The  ability  to  obtain  such 
a  covariance  must  be  assessed  for  the  specific  observation 
geometry  and  reentry  trajectory  being  estimated,  consistent 
with  the  linearization  assumptions  check. 

f.  Performance  is  improved  by  simultaneously  processing  ou^erva- 
tions  from  two  observation  satellites,  both  from  improved 
observability  and  higher  data  content,  unless  this  forces  rapid 
collapse  in  the  magnitude  of  the  state  covariance  matrix 
entries.  A  more  appropriate  observation  should  also  include 
range,  as  well  as  angular,  data. 

U.  With  mismatch  between  the  deterministic  estimator  dynamics  and  the 

truth  model  derived  data: 

a.  A  single  observation  per  update  epoch  is  recommended  to  insure 
local  model  validity. 

b.  At  least  a  minimum  number  of  observations  after  a  discrete 
change  in  state  variables  (i.e.,  a  step  change  in  3)  is 
necessary  to  insure  consistent  state  solutions. 

c.  In  its  present  infinite  memory,  deterministic  dynamics 
formulation,  significant  bias  exists  in  the  estimator 
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solutions..  The  next  chapter  will  address  this 
fundamental  limitation  of  the  estimator  for  application  to 
true  reentry  trajectories. 


Ill 


Chapter  III  -  Model  Compensation 
A .  Model  Compensation  Methods 

As  the  previous  numerical  results  indicated,  the  fundamental 
limitation  of  the  infinite  memory  estimator  formulation  with 
deterministic  dynamics  is  the  biased  estimator  solutions  which  occur 
when  processing  true  reentry  dynamics.  With  i)  exact  dynamics, 

ii)  an  upper  limit  on  the  time  span  of  valid  linearization,  and 

iii)  a  lower  limit  on  observation  accuracy  (greater  than  or  equal  to 
10~5  radians),  acceptable  estimator  performance  is  available  in  terms 
of  bias  and  RSS/(ONE  SIGMA)  ratio.  This  was  shown  in  the  Monte  Carlo 
results  of  Figures  6-13. 

If  significant  bias  had  existed  with  the  exact  dynamics, 
utilization  of  methods  which  address  the  character  of  neglected  higher 
order  terms  of  the  Taylor's  series  linearization  would  be  required. 
This  would  entail  examination  of  approximate  nonlinear  estimation 
methods  which  retain  selected  higher  order  terms  of  the  Taylor's 
series  expansions  for  the  reentry  application  (19,23,24).  Also,  one 
could  consider  extensions  of  a  statistical  linearization  approach, 
such  as  discussed  by  Gelb  (26)  for  the  scalar  dynamics  case,  f(x). 
These  methods  attempt  to  characterize  the  functional  form  of  the 
neglected  nonlinearities  in  a  polynomial  expansion  of  the  state  vari¬ 
ables.  Inherent  within  the  application,  however,  is  the  assumption 
that  a  valid  first  order  functional  description  of  the  true  dynamics 
exists  within  the  estimator  model,  f(x).  In  actual  application,  the 
variations  in  global  atmosphere  changes  and  reentry  vehicle 


fragmentation  effects  present  a  nearly  intractable  problem  for  mathe¬ 
matical  descriptions  of  all  reentries  within  a  deterministic  dynamics 
set. 

The  first  step  in  improving  the  linearized  estimator  performance 
concerns  the  impact  of  the  uncertainties  within  B  and  p.  One  must 
consider  model  compensation  methods  which  address  the  systematic  error 
resulting  from  the  inexact  and  uncertain  trajectory  dynamics.  The 
addition  of  pseudo-noise  to  the  state  dynamics  is  a  potential  means  to 
compensate  for  dynamic  uncertainties.  This  could  take  the  form  of  an 
additive  noise  term  to  the  state  dynamics  equation,  such  as: 

x(t)  =  f(x(t),t)  +  G(t)q(t)  (77) 

where : 

G(t)  is  a  time  dependent  coefficient  matrix, 

q  (t)  is  a  zero-mean  random  noise  term. 

In  the  linearized  application,  the  dynamics  equation  would  be 
reduced  to  the  variational  form: 

6x(t)  =  F(x  (t),t)  6x(t)  +  G(t)q  (t)  (78) 

o 

Incorporation  of  the  dynamic  noise  term  in  this  manner,  however, 
complicates  the  state  update  expression  (Equation  51).  It  also  requires 
an  extensive  tuning  process,  or  an  adaptive  technique,  to  modify  the 
coefficients  of  the  matrix  G(t)  to  incorporate  the  time  dependent 
effects  of  the  dynamic  uncertainties  properly  for  the  general  decay 
trajectory. 

An  alternative  to  this  formulation  is  discussed  by  Day  &  Schieb 
(27)  lor  application  to  the  differential  corrector.  This  aLso  involves 
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an  additional  noise  covariance  matrix  to  the  state  covariance  matrix 
during  the  propagation  phase.  In  this  approach,  the  state  covariance 
matrix  is  propagated  to  the  next  epoch  via: 


m+1  ,m 


=  ^ (tm+l ’ tm) 


[S  +  Q(t  ) J  ^(t  ,t  ) 
m  ,m  m  m+l  m 


(79) 


The  determination  of  the  Q(cm)  also  requires  extensive  simulation 
analysis  to  determine  the  noise  covariance  magnitudes  applicable  to  the 
widely  varying  dynamic  uncertainties  of  a  given  reentry  trajectory. 

This  approach  also  propagates  the  state  covariance  matrix  forward  based 
on  a  different  dynamics  model  than  the  state  trajectory.  Recall  that 
the  state  dynamics  propagate  according  to  Equation  2,  without  an 
additive  noise  term.  As  noted  previously,  Pon  (5)  has  pointed  out  that 
tuning  of  an  additive  pseudo-noise  matrix  for  orbiting  satellites  has 
proven  successful  only  when  significant  amounts  of  empirical  trajectory 
observations  are  available  to  aid  in  the  tuning  process.  This  pseudo¬ 
noise  matrix  was  not  transferable  to  other  satellites.  With  the  more 
uncertain  dynamics  and  the  short  trajectory  arcs  of  reentry, 
sufficient  empirical  data  is  not  available  for  each  decaying  satellite 
of  interest.  As  such,  a  more  systematic  approach  which  addresses  the 
time  dependent  character  of  the  dynamic  uncertainties  is  required. 

Statistical  adaptation  methods  offer  a  potential  means  to  address 
this  highly  time  dependent  nature  of  the  dynamic  uncertainties  for  a 
specific  decay  trajectory.  These  could  be  employed  in  an  approach 
similar  to  Jazwinski  (24)  to  allow  adaptive  estimation  of  noise 
covariance  terras,  or  in  the  more  general  sense  of  Maybeck  (25)  for 
determination  of  selected  uncertain  state  variables,  dynamic  noise 
matrix,  G(t),  or  the  dynamic  noise  covariance,  Q(t  ). 
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Unfortunately,  the  nature  of  the  current  application  precludes 
proper  incorporation  of  these  statistical  methods.  Statistical  adapta¬ 
tion  often  examines  the  character  of  non-zero  mean  residuals  or  an 
unbounded  growth  in  the  residuals  as  a  measure  of  when,  and  how  much, 
to  tune  the  noise  coefficients  (24). 

Alternately,  the  statistical  adaptation  methods  may  assume  a 
slowly  varying  state  variable  as  a  constant  over  n  observations.  Its 
value  is  then  determined  from  a  maximum  likelihood  formulation  for 
estimation  of  the  uncertain  state  variable  (25) ,  based  on  a  finite 
memory  of  measurements. 

In  either  case,  a  sufficiently  large  number  of  observations  are 
required  to  establish  the  statistical  validity  of  the  adaptation.  Due 
to  the  i)  rapidly  varying  dynamics  (particularly  b)  and,  ii)  the  need 
to  reduce  the  estimator  to  using  small  numbers  of  observations  at  an 
update  for  linearization  requirements,  these  methods  would  not  provide 
the  most  acceptable  adaptive  performance.  As  the  numerical  results  of 
this  chapter  show,  the  time  constant  of  the  variations  in  the  true 
dynamics  is  very  short  compared  to  the  10  second  data  rate  of  the 
current  observer.  Such  rapid  variations  in  true  dynamics  make  it 
difficult  to  assume  a  constant  value  over  the  required  n  observations 
for  statistical  analysis.  At  best  one  can  assume,  for  example,  6  is  a 
constant  between  successive  observation  points. 

Pon  (28)  has  suggested  a  technique  of  incorporating  an  additive 
noise  matrix  whose  terms  are  determined  from  a  single  orbit  estimation 
which  he  claims  will  avoid  the  need  for  a  priori  noise  covariance 
tuning.  Pon  alleges  his  technique  to  be  responsive  to  the  uncertain¬ 
ties  within  the  dynamics  between  successive  epochs  of  the  trajectory 


being  estimated.  At  a  given  epoch  t  ,  one  updates  the  state  via 
Equation  47  by  processing  the  next  k  observations  in  a  single  update. 
Simultaneously  at  the  next  epoch,  tm+^ >  one  obtains  a  second  estimate 
by  processing  the  same  k  observations,  now  at  times  prior  to  t  By 

propagation  of  the  t  estimate  forward  in  time  and  the  tm+^  estimate 
backward,  two  different  state  trajectory  solutions  exist  over  the  same 
time  span. 

By  trajectory  differencing  these  two  solutions  at  many  points 
along  this  time  span,  Pon  develops  a  matrix  of  second  order  statistics 
which  is  assumed  to  represent  the  dynamics  uncertainty  between  epochs. 
In  this  manner  he  obtains  a  diagonal  corruptive  noise  covariance  matrix 
for  insertion  into  Equation  79  to  propagate  forward  to  the  next  epoch, 
such  as: 

,  2 
x 

a  • 

x 


This  constructed  noise  covariance  matrix  suffers  from  several 
deficiencies.  It  would  require  a  large  number  of  numerical  computa¬ 
tions  to  develop  statistically  valid  o2  terms.  Also,  it  only 
compensates  the  diagonal  elements  of  the  state  covariance  matrix  at 

the  first  epoch,  t  .  In  this  "suboptimal"  sense,  it  neglects  the 
m 

cross-correlation  between  states  in  modifying  the  state  covariance 
matrix  to  represent  the  dynamic  uncertainties.  The  off-diagonal  terms 


are  modified  only  during  propagation  between  the  epochs.  A  more 


satisfying  technique  would  apply  a  full  pseudo-noise  compensation,  at 

the  update  epoch  and  during  propagation. 

Pon  (5)  shows  that  the  propagated  state  covariance  expression 

of  Equation  79  is  equivalent  to  premultiplying  by  an  upper  triangular 

form  matrix,  DCt^),  and  postmultiplying  by  a  lower  triangular  form 
T 

matrix,  D(t  )  ,  after  simple  state  covariance  propagation,  but  prior  to 
m 

update  as: 


S 

m,m- 


1 


=  <+> ( t  ,t  .)  [S  +  Q(t  )]  <f(t  ,t  )' 

m  m  1  m-l,ra-l  m-1  m  m-1 


=  D(t  )  * 4> ( t  >  t  )  S  <J>(t  ,  t  )T-D(t  )T 

m  m  m-i  m-i,m-l  m  m-1  m 


(81) 


This  technique  is  often  referred  to  as  a  "deweighting"  method.  Its 

effect  is  to  deweight  the  influence  of  the  output  of  the  dynamics 

model  when  S  .  is  used  in  the  next  epoch  update  (Equation  29) .  The 
m,m-l 

use  of  a  deweighting  matrix,  D(t  ),  is  an  extension  of  the  earlier 

m 

work  of  Fagin  (29)  ,  Sorenson  &  Sacks  (3) ,  and  Morrison  (2) ,  where  a 

scalar  deweighting  is  used.  With  a  scalar  deweighting,  the  influence 

of  the  old  observations  is  exponentially  deweighted  in  the  estimator 

solutions.  (Section  B  will  illustrate  this  effect  in  more  detail.) 

The  use  of  the  D(t  )  matrix  attempts  to  "deweight”  selected  states  by 
m 

differing  amounts  in  the  a  priori  state  covariance  matrix. 

In  some  applications,  the  D(t  )  matrix  has  been  applied  to 

m 

orbit  determinations  as  a  constant  deweighting  matrix  (5).  Similar 

to  tuning  the  noise  coefficients  within  Q(t  ),  the  diagonal  elements 

m 

of  D(t  )  are  developed  by  simulation  analysis, 
m 


D(t  ) 
m 


\ 


(82) 
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This  tuning  process  suffers  the  same  drawbacks  as  the  tuning  of 

the  noise  coefficients  in  Q(t  ).  Pon  (5)  states  that  translating  the 

in 

D(t^)  deweighting  matrix  from  one  satellite  to  another  proved 
unsatisfactory  for  orbital  applications.  Its  extension  to  the  more 
highly  variable  general  reentry  case  would  be  at  least  as  difficult, 
if  not  more  so.  Also,  acceptable  estimator  performance  was  achieved 
only  when  emperical  data  from  many  orbital  revolutions  were  available 
to  aid  in  the  tuning  process.  These  conditions  simply  do  not  exist 
for  the  reentry  application. 

Use  of  a  "suboptimal"  diagonal  deweighting  matrix,  however, 

provides  the  genesis  for  a  potentially  acceptable  adaptive  means  to 

determine  the  terms  of  a  deweighting  matrix,  D(t  ).  One  has  a  means  to 

m 

identify  the  potential  onset  of  bias  in  the  estimator  solution.  The 
standard  deviations  from  the  a  priori  state  covariance  matrix  reflect 
the  uncertainty  in  the  state  dynamics  and  previous  observation  history 
when  incorporated  into  the  state  update  Equations  29  and  47.  One 
assumes  local  dynamic  model  validity  over  the  time  span  between  the 
epoch  and  the  observation (s)  being  processed  for  an  update.  The 
magnitude  of  the  total  state  correction,  6x^  (from  Equation  47),  on 
each  state  can  be  compared  to  its  respective  standard  deviation  from 
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the  a  priori  covariance  matrix.  If  the  state  correction  is  greater  in 

magnitude,  one  can  consider  this  as  a  first  indication  in  the  growth  in 

bias  of  the  estimator  solution  which  is  inconsistent  with  the 

uncertainty  defined  by  the  state  covariance  matrix.  This  yields  a 

measure  to  determine  the  onset  of  a  bias  in  the  solution,  and  the 

direction  of  bias  development  along  a  particular  state  direction.  The 

advantage  of  this  measure  is  that  it  will  indicate  the  onset  of  a  bias 

in  the  solution  -  independent  of  its  cause.  Whether  the  bias  results 

due  to  dynamic  mismatch,  poor  observation  geometry,  or  excessive  noise 

on  the  observation  being  processed,  it  may  appear  as  an  excessively 

large  6x(t  ) .  in  one  or  more  state  variables  relative  to  their  a  priori 
m  1 

uncertainties. 

This  technique  was  chosen  for  its  simplicity  in  the  circumstances 
where  the  single  observation  updates  of  the  trajectory  are  required. 
Because  of  the  very  short  time  span  of  the  dynamics  model  validity,  an 
insufficient  number  of  observations  exist  at  the  epoch  update  point  for 
valid  statistical  testing  of  the  residuals.  A  valid  statistical 
determination  of  the  residual  mean  and  covariance  cannot  be  made  with 
a  single  sample.  The  comparison  of  the  6x^  estimate  to  its  a  priori 
standard  deviation  offers  no  improvement  in  this  statistical  quantifi¬ 
cation  since  the  same  information  available  in  the  single  observation 
is  used  to  develop  the  state  estimate.  It  does,  however,  provide  a 
mechanism  for  selection  of  the  scalar  deweighting  which  is  easily 
incorporated  into  the  estimator  structure. 

In  essence,  since  there  are  eight  states,  eight  individual 


conditions  exist  from  which  to  develop  a  model  compensation 
methodology.  Clearly,  the  use  of  these  eight  conditions  alone  do  not 


allow  one  to  determine  a  fully  populated  additive  noise  matrix  from  a 

single  trajectory,  or  a  multiplicative  upper  and  lower  triangular 

deweighting  matrix.  They  do  provide  sufficient  information  to  directly 

modify  the  eight  terms  of  a  diagonal  deweighting  matrix.  Therefore, 

initial  attempts  were  made  to  develop  the  diagonal  elements  of  the 

D(t  )  deweighting  matrix  (Equation  82). 
m 


Variations  on  the  implementation  of  this  deweighting  approach 
were  made  to  examine  the  potential  of  deweighting  only  selected 
variables  (i.e. ,  individual  states)  of  the  a  priori  covariance.  These 
included : 

1)  Deweight  only  8p  and  Q  (x_,x0)  by  a  constant  magnitude  at 
o  7  o 

each  epoch  throughout  the  trajectory  since  these  states 
contain  the  terms  which  most  directly  represent  the  dynamic 
uncertainties.  In  this  manner  the  deweighting  matrix  has 
the  form: 


D  ( t  ) 
m 


(83) 


2)  Similarly,  if  the  magnitude  of  the  total  state  correction  at 
any  given  epoch  is  greater  than  its  a  priori  standard 
deviation,  iteratively  select  the  affected  to  be  other 
than  unity,  leaving  the  other  d^  equal  to  one.  An  addi¬ 
tional  attempt  was  also  made  to  adaptively  deweight  only 
the  x?  and  xfl  terms  for  8p  and  Q,  leaving  positions  and 
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velocities  non-deweighted  (e.g.,  d^  =  1,  i  =  1  -  6). 

The  D(t^)  deweighting  matrix  is  used  to  modify  the  a  priori 

state  covariance  matrix,  S  .  from  Equation  81.  A  state  estimate  is 

m,m-l 

obtained  by  use  of  Equations  29  and  47.  This  process  repeats,  until: 

6x(t  ).  <  /~S  "  i  =  1,8  (84) 

mi-  m ,m-l . . 

li 

Both  attempts  proved  unsuccessful.  The  primary  causes  of  these 
methods  failing  to  improve  the  biased  estimator  performance  are  the 
following : 

1)  This  is  clearly  a  suboptimal  approach  to  deweighting.  The 

diagonal  deweighting  matrix,  D(t  ) ,  is  not  equivalent  to 

m 

applying  a  full  pseudo-noise  matrix  to  the  state  covariance 

matrix.  A  triangular  form  of  D(t  )  would  be  required  for  a 

m 

full  pseudo-noise  compensation. 

2)  The  artificial  deweighting  of  selected  states,  and  not 
others,  improves  the  BIAS/ (ONE  SIGMA)  ratio  on  the 
deweighted  states  at  the  expense  of  aggravating  the  bias  on 
the  other  states. 

3)  The  deweighting  on  only  the  x,  and  xQ  states  by  a  constant 
magnitude  throughout  reentry  was  not  responsive  to  the  time 
dependent  nature  of  the  dynamic  uncertainties. 

As  a  consequence,  the  mismatch  between  state  corrections,  <$x^ , 
and  their  a  priori  standard  deviations  was  used  to  determine  a  simple 
scalar  deweighting  of  the  a  priori  covariance  matrix.  This  scalar 
deweighting  can  be  considered  as  a  special  case  of  the  diagonal 
deweighting  where  all  the  d^  terms  are  equal.  This  scalar  deweighting 
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offered  superior  performance  than  the  attempts  to  deweight  the  individ¬ 
ual  states  by  different  amounts  for  a  number  of  reasons. 


1)  It  is  strictly  derivable  as  an  exponential  decay  of  the 
filter  memory.  (2,3,29,30) 

2)  As  the  numerical  results  presented  later  will  demonstrate, 
it  avoids  the  requirement  for  extensive  simulation 
analysis.  The  comparison  of  the  6x^  magnitudes  with  their 
a  priori  standard  deviations  successfully  determines  the 
need  for  deweighting.  With  proper  selection  of  the  scalar 
deweighting  constant,  acceptable  estimator  performance  is 
evident  which  responds  to  the  time  dependent  nature  of  the 
dynamics  mismatch.  As  the  Monte  Carlo  results  show,  the 
deweighted  state  covariance  matrix  provides  a  satisfactory 
measure  of  the  RSS  uncertainty  in  the  estimator  position 
and  velocity  solutions,  and  the  bias  in  the  state  estimate 
is  less  than  the  a  priori  standard  deviation  from  the 
deweighted  state  covariance  matrix. 

The  Monte  Carlo  results  on  the  basic  estimator  with  mismatched 
dynamics  (Figures  14-17  of  the  previous  chapter)  substantiate 
acceptable  estimator  performance  in  the  early  trajectory  phases,  in 
terms  of  bias  and  RSS/(0NE  SIGMA)  ratio  relative  to  the  estimator 
computed  standard  deviations.  This  results  since  the  deterministic 
estimator  dynamics  model  can  locally  represent  the  truth  model 
dynamics  in  this  region.  The  truth  model  dynamics  are  adequately 
approximated  by  the  exponential  atmosphere  and  constant  gp^  between 
epochs.  Secondly,  a  small  number  of  observations  have  been 
processed.  The  terms  of  the  state  covariance  matrix  have  not  reached 
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such  a  small  magnitude  that  the  output  of  the  deterministic  model 
dynamics  becomes  overly  weighted  relative  to  the  information  matrix  of 
the  new  observations. 

The  more  dramatic  dynamics  variations  occur  later  in  the 
trajectory  data  span.  In  this  region  the  magnitude  of  the  terms  in 
the  state  covariance  matrix  continue  to  decrease  in  the  velocity  terms 
and  begin  to  also  decrease  in  the  position  terms.  The  output  of  the 
deterministic  estimator  model  becomes  increasingly  weighted,  thereby 
causing  the  estimator  to  begin  to  ignore  the  complete  information 
content  in  the  later  observations. 

The  fading  memory  estimator  alleviates  much  of  this  phenomena. 
It  allows  greater  flexibility  of  the  estimator  to  accommodate  the  new 
observations  by  increasing  the  magnitude  of  terms  in  the  a  priori 
state  covariance  matrix  at  each  new  epoch  update  point.  The  effects 
of  this  scalar  deweighting  technique  are  shown  in  Section  D.2.  of  this 
chapter  with  the  presentation  of  the  results  of  the  Monte  Carlo 
analysis. 

An  alternate  form  of  dynamic  compensation  was  also  considered 

which  would  combine  the  use  of  an  additional  pseudo-noise  matrix  and 

the  scalar  deweighting  or  "tuning"  parameter,  y,  as  Q(t)  = 

Q  (t)  +  y  A  Q(t).  One  may  adaptively  determine  the  magnitude  of  the 
o 

scalar  tuning  parameter,  perhaps  from  the  techniques  discussed  above. 
The  problem  still  remains  to  construct  the  pseudo-noise  matrix 
properly  at  the  initial  reentry  epoch  point.  In  essence,  the  proposed 
formulation  is  analogous  to  incorporating  the  Q(t)  =  ( t )  +  y  A  Q(t) 

into  the  estimator  structure.  In  the  absence  of  more  adequate  infor¬ 
mation  to  develop  the  initial  Qq  and  AQ  terms,  the  estimator-computed. 
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or  a  priori,  state  covariance  matrix  provides  a  description  of  the 


uncertainties  which  is  modified  at  succeeding  epochs  by  the  adaptive 
selection  of  the  scalar  parameter,  y.  Section  B.  ,  below,  will  show 
the  details  of  the  fading  memory  estimator  development. 

B .  The  Fading  Memory  Differential  Corrector 

Based  on  an  exponential  decay  of  the  observation  weighting, 
several  authors  have  extended  the  early  work  of  Fagin  (29)  for  the  age 
weighting  of  scalar  observations.  Tarn  and  ZaborcVv  (30)  report  an 
extension  of  age  weighting  to  the  vector  observation  case  in  1970. 
Sorenson  and  Sacks  (3)  develop  an  exponential  aging  technique  for  the 
observations  in  a  linear  Kalman  Filter  formulation.  They  also  point 
out  that  the  exponential  aging  applies  to  three  distinct  quantities 
within  the  filter  structure:  i)  the  initial,  a  priori  state 
covariance  matrix,  li)  the  observation  covariance,  (previous 

observations) ,  and  iii)  the  state  covariance  matrix  prior  to  update 
at  each  new  observation.  Morrison  (2)  shows  a  derivation  for  the 
exponential  aging  from  the  initial  weighted  least-squares  cost 
function  for  both  a  Bayes  Filter  and  Kalman  Filter  for  linear  applica¬ 
tions  . 

The  following  derivation  (Equations  85-103)  will  directly  follow 
Morrison's  (2)  for  a  linear,  fading  memory  estimation  technique  which 
is  now  applied  to  the  differential  -  corrector  estimation  of  the  non¬ 
linear  reentry  problem.  The  initial  presentation  will  concentrate  on 
the  batch  processor.  This  formulation  will  then  be  modified  to  yield 
a  recursive  formulation  of  the  fading  memory  estimator  which  is 
applicable  to  the  uncertain  dynamics  of  the  reentry  trajectory.  The 
differences  in  this  development  from  Morrison  (2)  are  only  due  to  the 
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iterative  updates  of  the  linearized  differential  corrector  and  the 
development  of  an  alternate  expression  for  the  ”non-deweighted"  state 
covariance  matrix  expression,  which  may  be  obtained  from  the  deweighted 
covariance  matrix  expression. 

One  may  define  a  modified  cost  function  as  follows: 


J  "  t\„)  -  Tfa>  55<VlT  !“(„)  -  T(»)  ftV'MV,)!"  <85) 


where  the  scalar  multiple,  y"1  =  [y(t,  ,t  )]~1  takes  on  values: 

k  n 


0  <  1/v  <  1 


(86) 


based  upen  the  age  of  the  observations.  In  the  batch  formulation  k 

takes  on  the  values  k  =  l,...,n-l,  and  v(t  .t  )  =  1,  at  the  time  of  the 

n  n 

latest  observation.  Recall  that  t  is  the  epoch  time,  v,  .  is  the 

o  K  (n) 

vector  of  observation  residuals  (Equation  18),  and  the  matrix 

incorporates  the  observation  geometry  and  the  state  transition  matrix 
between  the  epoch  and  the  observations  (Equation  16) . 

In  matrix-vector  formulation,  the  weighting  matrix,  R^,  on  the 
observations  may  be  redefined  to  reflect  the  accumulation  of  the 
deweighting  scalar  in  a  slightly  altered  f 


R  •  1 
(  n 


R  ,  Y  ( t  >  t  ) 
n-1  n  n-1 


R(n) 


R1  Y(tn*tl)  ' 


(87) 
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where  n  is  the  observation  at  the  current  epoch  being  updated,  and  the 
Y  =  y(t  ,t,  )  reflects  increasing  uncertainty  on  the  older  observations. 
In  this  manner,  an  exponential  memory  is  established  (2)  when  the 
[^R(n)  (y)J  inverse  is  incorporated  into  the  estimator  structure  to 
weight  the  observations. 

The  cost  function  may  now  be  written  in  the  matrix-vector  form 
where  the  argument  on  y  has  been  neglected: 

J  ‘  -  TM  -<‘o»T[V>  W]"1  -  T(n)  5;<to)!  t88) 

As  before  since  (y)J-1  Is  positive  definite,  the  necessary  and 

sufficient  conditions  for  minimizing  J  are: 


3  J 

3<5x(t  ) 
o 


(89) 


Making  use  of  the  matrix  identity, 
block  diagonal  weighting  matrix,  |r 
minimization  of  J  results  in: 


(AB)T  = 


T  T 

bV 


and  the  symmetry  of  the 


the 


(n) 


T(n)6x(to) 


=  0 


(90) 


Solving  for  6x(to)  provides  the  batch  processing  update 


expression : 


6x(to)  =  (T 


(91) 


where  now  the  "deweighted"  state  covariance  matrix  after  update  is 
obtained,  assuming  a  positive  definite  information  matrix. 
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(T(n>  [R(»>  W]  “  T<»> : 


p„.„  ■  «<„)  [r<„>  *<»>>•■  <92> 

Note  that  in  this  batch  formulation,  the  covariance  matrix,  P  ,  from 

n  ,n 

Equation  92  assumes  no  existing  a  priori  covariance  matrix  information, 
but  is  simply  derived  from  the  estimator  processing  of  the  observations. 
The  non-deweighted  state  covariance  matrix  can  also  be  derived.  By 
definition,  for  zero  mean  Sx(tQ): 


S  =  E(6x(t  )  5x(t  )  ) 
n,n  o  o 


where  6x(tQ)  is  obtained  from  Equation  91. 

Substitute  for  6x(t  ) : 

o 


S  =  E{ 
n,n 


[(T(n,T[R(„)M]“  T<„>>"  T(„)T[R(,)  <''>]“  P(n,l  <«> 


[V)T(%>W]“ 


T(n)rl  T(n)T[R<n)W]‘1  v(n),T' 


Taking  the  transpose  and  simplifying  yields: 


Sn.„  ’  <T<„)  [R<n>«]'‘  V‘  T(„)  [R(„)«]'1  «*<„>*(„> 


E(v,  vV,  N  )  (95) 


[R<n)<»]'1  WT(n)  [RC„>«]“ 


where 


[%)«]“ 


and  (T 


<n)  [R(„>W]“  T(n) 


T.  .)“  are  positive  definite 


matrices . 


Using  the  definition  for 


[R(»)w]‘‘ 


and  the  definition 


R.  .  =  E(v,  ,v.  .  ),  the  middle  term  of  Equation  96  becomes: 
fnj  (n)  (nj 

[•wWr-Mf'wW)'1  ' 

■  <R<„)<^>r‘ 
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Therefore,  the  non-deweighted  covariance  can  now  be  related  to  the 
deweighted  covariance  by: 


S 

n  ,n 


(Y) 


=  P 

n,n 


(Y) 


(97) 


Remember,  the  S  matrix  is  a  valid  minimum  variance  expression 
n  ,n 

under  the  assumption  that  <5x(tQ)  is  zero  mean  -  or  that  the  estimate  of 
the  state  yields  an  unbiased  solution.  This  is  an  important  considera¬ 
tion  which  must  be  examined  in  the  later  numerical  simulations.  If  the 
5x(tQ)  is  not  an  unbiased  solution,  it  may  be  imprudent  to  use  the 

second  central  moments  S  statistics  for  a  measure  of  the  error  in 

n,n 

the  estimator  solution. 

Now  consider  a  recursive  formulation  where  each  epoch  is 

designated  by  the  index  m.  The  observations  are  designated  by  the  index 

n,  where  n=l,...,k.  One  can  partition  the  T^  and  ^(n) (y)  matrices  to 

isolate  the  latest  observation  at  time  t  and  establish  the  update 

mn 

epoch  t  between  the  t  and  t  observation  times.  Recall  that  H  is 
m  m  ,  m  n 

n-1  n 

defined  in  Equation  14. 

t  \ 


T 


(n) 


H  *<t  ,t  ) 
n  mm 
_ n _ 

T,  n  <Kt  ,t  ) 
(.n-1)  m  ,  m 
n-l 


(98) 


and 


R 

n 


W'1 


< 


R,  . , [ y ( t  ,t 
(n-1)  n  n- 


,)] 


-i 


(99) 


V  i  > 

By  substitution,  and  after  simplification,  the  inverse  of  the  updated 

deweighted  covariance  at  epoch  time,  t  ,  becomes  (2) 

m 
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P«,»'1  ■  T<n)T[R(„)(rt]'‘  T(n) 

'  TnT  V*  T„  +  lh  ♦(Wl> 

•  T,  ..  4>(t  ,t  .)T 
(n-1)  m  m-1 

or 


(100) 


P  _1  =  (1/y  p  ":  +  t  T  R  _1  T  )  (101) 

m,m  m,ra-l  n  n  n 


where  the  a  priori  deweighted  state  covariance  is  obtained  by  propagat¬ 
ing  forward  from  a  previous  epoch  at  time  t^_^.  This  epoch  contains 
the  state  estimate  for  the  observations  numbered  1  through  n-1. 

In  recursive  form,  the  state  update  on  the  next  (non-deweighted) 
observation(s)  is  now: 


«x(t  )  =  (1/y  Pm  m  r1  +  T  T  R  -1  T  )“ 1  T 1  R  " 1  V 

m  m,m-l  n  n  n  n  n  n 


(102) 


where: 

Y  is  now  a  one-step,  deweighting  scalar 

P  ,  =  a  priori  deweighted  covariance 
in, m-1 

As  in  the  infinite  memory  case.  Equation  102  is  applied  iteratively  at 

each  epoch  until  the  convergence  criteria  are  satisfied.  The 
X  T 

T  R  _1  T  and  T  R  -1  v  matrices  are  evaluated  from  the  reference 
n  n  n  n  n  n 

trajectory  on  the  final  iteration,  with  R^  the  covariance  of  the 
observations  at  time  t 

m 

n 

The  updated,  deweighted  state  covariance  matrix,  P  ,  is  not 

m,m 

the  non-deweighted  covariance  matrix,  S  .  But  it  does  represent  the 

ni  >  rn 

model  uncertainty  and  the  scalar  deweighting  of  the  old  observations  at 
each  epoch  point  and  is  so  incorporated  into  Equation  102.  This 
provides  the  recursive  differential  corrector  formulation  which  allows 
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selection  of  y  at  each  epoch  to  reflect  the  uncertainty  in  the  previous 

history  of  the  dynamics  model. 

As  the  numerical  results  of  Section  D.  show,  the  use  of  the 

P  covariance  matrix  yields  a  conservative  estimate  of  the  errors  in 
m  ,m 

the  state  solution.  The  use  of  the  non-deweighted  covariance,  S  ,  is 

m  ,m 

not  necessary  for  on-line  use  of  the  estimator  for  this  application. 

It  may,  however,  be  implemented  into  the  estimator  structure.  Should 

improvements  in  the  dynamics  model  become  available  such  that  zero  mean 

state  estimates  can  be  obtained,  the  S  matrix  would  provide  a  minimum 

m,m 

variance  covariance  matrix  from  the  estimator  computations. 

Let  the  development  now  depart  from  following  that  of  Morrison 

(2) ,  where  he  shows  an  alternate  relationship  between  S  and  P  in 

m,m  m,m 

Equation  97.  Doing  so,  he  develops  a  recursive  expression  for  S 

m  ,m 

and  S^  which  will  not  be  repeated  here.  Morrison's  expression  was 
developed  for  a  linear  Bayes  and  Kalman  Filter  application.  In  its 
differential  corrector  form,  the  update  expression  for  the  non- 
deweighted  covariance  at  each  new  observation  would  be  structured  as 
follows : 


S  =  (I  -  P  TTR_1T)S  .  (I  -  P  t  T  r  '  1 
m,m  m,m  n  n  n  m,m-l  m,m  n  n 


+  P  T  T  R  "  XT 
m,m  n  n  n 


P 

m,m 


(103) 


Due  to  the  highly  ill-conditioned  nature  of  the  angle  only 

X 

derived  single  observation  information  matrix,  T  R  _1  T  ,  and  the 

n  n  n 

lack  of  positive  definiteness  for  single  observation  updates  discussed 

earlier,  this  formulation  was  found  to  be  numerically  unstable.  The 

T 

subtraction  of  P  T  R  " 1  T  from  the  identity  matrix,  I,  involves 
m,m  n  n  n 
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terms  which  are  very  close  to  zero  or  to  one.  The  subtraction  of  these 

nearly  identical  numbers  results  in  an  accumulation  of  error  when 

implemented  on  the  finite  word-length  digital  computer.  These  errors 

are  further  aggravated  when  ^  is  propagated  forward  to  new  epochs  by 

use  of  the  state  transition  matrix,  <f>(t,t  ).  This  instability  results 

m 

from  the  limited  observability  of  the  single  observer  data  used  to 

T 

formulate  the  information  matrix,  T  R  - 1  T  .  Recall  that  Section  D.3. 

n  n  n 

in  Chapter  II  described  the  numerical  difficulties  associated  with  use 

of  the  single  observer  angular  data. 

To  correct  for  this  implementation  problem,  consider  the 
T 

(I  -  P  T  R  -  1  T  )  expression.  From  Equation  101,  one  may  solve 
m,m  n  n  n 

for: 


T  R  -1  T  =  P  -1-  1/yP  -1 
n  n  n  m,m  m,m-l 


(104) 


Substituting  into  the  [I-P  T  R-1T]  expression  yields: 

m,m  n  n  n 


[I-P  (P 


- 1 


1/Y  p  .-1)]  =  [I  -  I  +  P  1/y  p 


m,m  m,m 


-  1 


‘m,tn-l 


m,ra  m,m-l 


]  (105) 


=  [p  i/y  p  -1: 

m  ,m  m ,m- 1 


Substituting  back  into  the  S  expression  of  Equation  103, 

rn ,  m 


s  =  [P  1/y  p  "M  s  ,  [I/y  p  -1  p  ] 

m,m  m,m  m,m-l  m,m-l  m,m-l  m,m 

+  P  ttr~1t  P 
m,m  n  n  n  m,m 


one  obtains: 

006) 


which  provides  a  more  stable  numerical  expression  which  allows  update  of 
the  non-deweighted  covariance  matrix  at  the  new  observation  epoch. 


The  illustration  below  will  show  the  relative  stability  characteristics 
of  both  Equations  103  and  106  for  the  total  position  standard  deviations 
from  a  typical  estimator  run. 
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With  the  above  derivations  for  the  state  update  and  deweighted 
covariance,  a  complete  listing  of  the  recursive  form  of  the  fading 
memory  differential  corrector  may  be  shown.  Each  epoch  point  will  be 
designated  with  the  index  m.  In  general,  multiple  observation  updates 
may  be  made  at  an  epoch,  therefore  the  observations  are  designated  by 
the  index  n,  where  n=l,...,k.  In  application  to  the  reentry 
trajectory  with  uncertain  dynamics,  n  often  equals  one,  a  recursive 
formulation  for  each  observation.  Recall  that  the  state  update 
equation  is  applied  iteratively  until  convergence  on  iteration  i,. 
Propagation  Between  Epochs: 

State :  via  integration  of  x(t)  =  f(x,t)  (107) 

Covariance : 

Deweighted : 


P 

m,m- 


1 


4>(t  ,t  .) 

m  m-1 


m-1 ,m-l 


4>(t  ,t  A 

m  ra-i 


(108) 


Update  at  New  Epoch: 


State:  x(t)=x  .  (t  )  +  E  <$x(t  ). 

- —  mm  m-lm  ...  mi 

1=1 


6x(t  )  .  =  (1/y  P  "1+TiR-iT)-iTiR  '  v 
m  i  m, m-1  n  n  n  n  n  n 


(109) 

(110) 


Covariance : 

Deweighted : 

P  =  (1/y  P  “x  +  T  T  R  -x  T  )-1  (111) 

m,m  m,m-l  n  n  n 

One  may  also  compute  the  non-deweighted  covariance  matrix  with  the  on¬ 
line  estimator  software,  if  desired: 
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Propagation : 


S 

m,m- 


1 


<f>(t  ,t  ) 
m  m-1 


S  .  .  4>(t  ,t  .) 

m-1, m-1  m  m-1 


Update : 


(112a) 


S 

m  ,m 


[p 

m,m 


1/y  P 


-  1  ■ 


m,m-l 


S  i [ 1/y  P  . 
m,m-l  m,m-l 


- 1 


P  ] 

m,m 


T 

+  P  T 
m,m  n 


R  -1 
n 


T  P 
n  m  ,m 


(112b) 


Application  of  the  estimator  is  similar  to  the  use  of  the 

infinite  memory  estimator  where  small  numbers  of  observations  must  be 

processed  at  each  epoch  update  point.  The  number  of  observations 

allowed  are  determined  from  the  linearity  check  of  Equation  44.  The 

fundamental  difference  is  to  select  the  scalar  deweighting  constant, 

Y,  such  that  sufficient  fading  occurs  to  allow  the  estimator  solution 

to  adjust  to  the  newer  observations.  Proper  selection  of  Y  allows  the 

estimator  to  deweight  the  output  of  the  dynamics  model  via 

(1/y  P  ,-1)  at  each  new  epoch, 
in  jin — JL 

Both  Morrison  (2)  and  Sorrenson  and  Sacks  (3)  suggest  selection 
of  the  scalar  constant  by  analysis  of  the  residuals  of  the  solution 
process.  As  noted  earlier,  they  also  point  out  that  observation 
residual  testing  is  most  valid  when  i)  the  functional  form  of  the 
dynamics  modal  is  fundamentally  sound,  and  ii)  sufficiently  large 
numbers  of  residuals  are  available  over  this  span  for  statistical 
analysis  of  the  residuals  to  be  valid.  These  conditions  do  not  exist 
for  this  reentry  application.  For  the  reentry  application,  local 
model  validity  is  assumed  between  epoch  and  the  observation  being 
processed  in  a  single  or  recursive  formulation. 
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If  the  5x. (t  )  <  /  P  ,  criterion  is  violated  at  any  given 
i  o  -  m,m-I^ 

epoch,  y  may  be  chosen  such  that  the  resulting  correction  to  any  given 
state  variable  is  within  the  one  sigma  level  of  the  a  priori  deweighted 
covariance  matrix.  This  enhances  the  consistency  in  the  state  estimate 
with  the  uncertainty  in  the  previous  estimator  solution,  and  provides 
an  ad  hoc  technique  to  determine  the  need  for,  and  magnitude  of,  the 
fading  memory  required.  It  also  allows  more  freedom  in  the  state 
update  by  processing  a  new  observation  with  a  lesser  weighting  on  the 
output  of  the  estimator  dynamics  model  or  history  of  previous 
observations.  This  adaptive  selection  of  a  fading  memory  parameter  is 
very  much  like  the  scalar  tuning  parameter  within  a  noise  matrix 
Q(t)  =  QQ(t)  +  y  A  Q(t)  should  the  initial  Qq  and  AQ  be  capable  of 
being  constructed. 

A  principal  consideration  of  this  adaptive  selection  of  a  fading 
memory  was  to  avoid  the  need  for  extensive  simulation  analysis  for 
tuning  an  additive  pseudo-noise  matrix.  The  numerical  results  which 
follow  later  will  illustrate  acceptable  estimator  performance  (bias 
and  RSS/(ONE  SIGMA)  ratio)  is  achieved  relative  to  the  deweighted 
state  covariance.  However,  there  does  exist  some  sensitivity  in  the 
selection  of  this  ad  hoc  scalar  tuning  parameter:  One  cannot  allow  a 
large  y  to  be  selected  early  in  the  trajectory  before  the  reentry  is 
committed,  due  to  the  stability  considerations  discussed  below. 
Fortunately,  the  dynamics  uncertainties  and  observation  accumulation 
which  contribute  to  biased  estimator  performance  generally  do  not 
require  a  fading  memory  in  these  very  early  trajectory  phases. 


C.  Stability  Considerations 
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A  rigorous  stability  analysis  lias  not  been  accomplished  for  the 
proposed  estimation  algorithm.  Potential  applications  of  the  fading 
memory  estimator  must  consider  this  limitation.  The  following  discus¬ 
sion  lends  some  insight  into  the  question  of  estimator  stability. 

The  question  of  estimator  stability  is  a  complex  one  for  the 
current  application.  The  stability  characteristics  are  dependent  on 
the  underlying  stability  of  the  reentry  dynamics  process.  Miller  (29), 
Morrison  (2),  and  Sorrenson  and  Sacks  (3)  all  discuss  the  stability 
aspects  of  their  respective  linear  estimator  formulations  which  use  a 
fading  memory.  The  three  discussions  center  around  the  linear 
character  of  the  estimator,  allowing  reasonably  straightforward 
applications  of  linearized  stability  theory. 

Miller  (29)  assumes  a  continuous,  stationary  system  with  linear, 
constant  coefficients  while  examining  the  stability  of  the  inverse 
covariance  propagation  via: 

d/dt(P(t)_1)  =  a  P(t)-1  +  FP(t)-1  +  P(t)"1  FT  (113) 

-  P(t)"1  HT  R'1  H  P(t)-1 

where: 

a  =  a  constant  fading  memory  parameter 
F  =  a  constant  dynamics  matrix 
H  =  a  constant  geometry  matrix 
R  =  a  constant  observation  covariance  matrix 
Also,  in  a  continuous  formulation,  Sorrenson  and  Sacks  (3) 
analyze  the  stability  of  a  linear  (uniformly  and  completely  controllable 
and  observable)  system.  They  relate  their  findings  to  the  asymptotic 

136 


stability  characteristics  of  the  infinite  memory  Kalman  Filter. 

Morrison  (2)  discusses  the  stability  aspects  of  a  discrete 
observation  formulation  for  a  fading  memory  estimator  by  examining: 


P  _1 
ra,m 


(1/y)  <J>(t  >  t  1)"T  P  .  t"1  <Kt  ,t  .r1  +  ht  R"1  h 

'  T  m  m-1  m-l,m-l  T  m  m-1 


(114) 


which  assumes: 

1)  a  constant  coefficient,  linear  model 

2)  a  constant  coefficient,  linear  observation  relationship 

3)  a  constant  input  covariance  matrix 

4)  a  constant  stepsize  between  observations 

The  basic  tenet  of  any  of  these  approaches  is  to  show  that  the 
state  covariance  matrix  does  not  grow  unbounded  in  an  asymptotic  sense; 
and  thereby,  drive  the  state  solution  unbounded.  In  the  current  reentry 
application  the  covariance  matrix  propagation  and  update  equations 
include  both  time  dependent  variables  and  a  nonlinear  combination  of 
state  variables.  A  linearized,  constant  coefficient;  linearized, 
periodic  coefficient;  or  perturbations  theory  stability  approach  is 
not  appropriate.  The  time  dependent  and  nonlinear  dynamics  are  the 
reason  that  these  techniques  are  not  appropriate. 

The  structure  of  the  current  reentry  estimator  yields  two  facets 
which  merit  investigation  regarding  the  stability  of  the  algorithm. 

One  is  the  state  update  expression  which  processes  an  observation  to 
modify  the  state  estimate  at  an  epoch.  It  is  desirable  to  examine 
the  nature  of  the  bound  on  the  resulting  state  estimate  as  a  function 
of  the  observation  data,  the  value  of  the  fading  memory  parameter,  and 
the  a  priori  state  covariance  matrix  being  used  to  update  the  state 
estimate.  The  other  stability  consideration  concerns  the  propagation 


of  Che  state  trajectory  forward  from  an  epoch  by  integration  of  the 
system  dynamics  equation.  Two  different  approaches  to  these  stability 
questions  were  considered.  Safonov  (32)  offers  a  functional  analysis 
approach  which  may  potentially  be  used  to  examine  the  numerical 
stability  of  the  system  during  the  state  update  phase.  The  propaga¬ 
tion  phase  may  be  better  suited  to  the  energy  considerations  of  the 
system  dynamics  using  an  asymptotic  stability  analysis  of  Liapunov  (33) . 

Safonov's  approach  (32)  suggests  a  topological  separation  of 
spaces  to  identify  regions  of  stable  and  unstable  behavior  relative  to 
the  reference  trajectory  state  solution.  The  state  update  equation 
and  observation  equation  may  be  cast  into  the  form  of  a  multivariable 
feedback  system.  Each  individual  iteration  of  the  estimator  at  an 
epoch  update  point  could  be  written  in  such  a  form.  Since  the 
linearized  equations  are  evaluated  along  the  reference  trajectory, 
constant  values  of  the  geometric  partials  matrix  and  the  state 
transition  matrix  are  available  for  examination  at  a  given  iteration. 

The  Safonov  approach  is  a  higher  level  of  abstraction  which 
includes  the  Liapunov  stability  approach  for  the  multivariable  feedback 
system.  However,  the  stability  characteristics  of  the  system  are 
dependent  on  the  adequacy  of  the  linearization,  the  accuracy  of  the 
estimator  dynamics  relative  to  the  true  dynamics,  and  require 
continuous  partial  derivatives  of  the  f(x)  term.  The  reentry  problem 
has  uncertain  true  dynamics,  does  not  provide  an  exact  linearization, 
and  does  not  have  continuous  partial  derivatives  for  the  true  dynamics 
due  to  vehicle  fragmentation,  although  the  estimator  dynamics  model 
does  have  continuous  partial  derivatives  of  f(x).  Substantial  research 
would  be  required  to  adapt  the  Safonov  approach  to  the  current 
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application.  Unfortunately,  the  examination  for  stability  would  also 
need  to  be  accomplished  for  each  iteration  at  each  update  epoch  due 
to  the  re-linearization  about  a  new  reference  trajectory  as  the 
estimator  converges. 

The  propagation  of  the  state  estimate  forward  from  an  epoch  may 
occur  over  regions  where  no  observations  are  available  and;  therefore, 
may  not  be  suitable  for  expressions  of  the  multivariable  feedback  form. 
The  energy  considerations  of  the  system  dynamics  may  potentially  be 
examined  for  their  asymptotic  stability  characteristics  using  a 
Liapunov  approach  (33) .  One  could  potentially  examine  the  asymptotic 
stability  characteristics  of  the  state  trajectory  solution  at  each 
epoch  update  point.  A  particular  selection  of  the  adaptively 
determined  fading  memory  parameter,  y,  in  combination  with  the 
asymptotic  stability  characteristics  of  x(t)  =  f(x(t),t)  could  be 
examined  at  each  epoch  relative  to  the  EC1  coordinate  system  origin, 
the  Earth  center. 

The  current  dynamic  system  does  not  consider  any  rotational 
energy  or  motion,  but  is  simply  a  point  mass  formulation  for  the  reentry 
trajectory  dynamics.  This  makes  a  formulation  similar  to 
Meirovitch  (33)  appealing.  The  system  kinetic  energy,  T,  potential 
energy,  V,  and  the  dissipative  drag  forces  were  considered  to 
formulate  a  system  Hamiltonian,  H',  for  use  as  a  possible  Liaponov 
function.  Meirovitch  (33)  shows  the  ability  to  construct  a  Rayleigh 
dissipative  function  which  is  a  quadratic  function  in  the  velocity 
terms  when  the  non-conservative  dissipative  forces  are  linearly 
dependent  only  on  the  velocity  terms: 
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f  =  -c  x  (115) 

XX 

The  Rayleigh  dissipative  function  takes  the  general  form: 

F  =  *5  H  c,.  q.  q,  (116) 

i  j  i]  1  J 

where 

c..  =  constant  coefficients 
11 

=  generalized  velocity  terms 

In  this  manner,  the  system  Hamiltonian,  H’  =  T  +  V,  and  the 
quadratic  form  of  the  total  time  derivative  of  H'  are  examined,  where: 

H'  =  -2F  (117) 

If  H'  is  positive  definite  and  H'  is  positive  semidef inite ,  the  system 
Hamiltonian  may  satisfy  one  of  the  basic  Liapunov  stability  theorems 
and  a  statement  on  its  asymptotic  stability  can  be  made.  Specifically, 
If  the  set  of  points  where  H'  =  0  contains  no  nontrivial  solution 
trajectory  for  the  dynamic  system  x(t)  =•  f(x,t),  the  solution  is 
asymptotically  stable  (33). 

Unfortunately,  the  dissipative  forces  in  the  current  estimator 
model  are  more  complex  in  the  velocity  terms  and  coupled  to  the 
position  terms  due  to  the  rotating  Earth  atmosphere  as: 

Ed  =  m  Vr  (x  +“>y) 

X 

fj  =  -h  m  Bp  VR  (y  -u>x)  (118) 

y 

f  d  =  m  Bp  VR  z 
z 

where:  t.i  =  Earth  rotation  rate,  and  VR  =  [(x+my)2  +  (y-mx)2  +  z2]^. 


-  Ti  1 1  !■ 
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While  the  system  Hamiltonian,  H'  =  T  +  V,  is  positive  definite 


in  the  velocity  terms,  the  total  time  derivative,  H',  is  not  quadratic 
in  the  velocity  terms.  It  must  therefore  be  examined  at  every  point 
along  the  reentry  trajectory  to  gain  any  insight  into  the  system 
stability  characteristics.  To  accomplish  this,  one  must  integrate  the 
system  dynamics  equation  for  each  epoch  solution.  There  is  no  apparent 
quadratic  form  whose  sign-definiteness  could  be  examined  to  identify 
the  stability  characteristics  of  the  system.  Early  during  reentry,  for 
a  given  value  of  y  and  observation  geometry,  the  reentry  dynamics  may 
indeed  become  unstable.  The  resulting  state  trajectory  solution  will 
skip  the  Earth  atmosphere. 

From  strict  physical  considerations,  one  can  restrict  the 
estimator  application  in  these  upper  altitude  regions.  The  magnitude 
of  the  total  velocity  at  an  epoch  can  be  examined  relative  to  the 
escape  velocity  at  that  particular  altitude  (34): 


V  = 
e 


G (M+m) 


(119) 


where:  G  =  the  universal  gravitational  constant 

M  =  the  mass  of  the  Earth 
m  =  the  mass  of  the  satellite 

R  =  the  magnitude  of  the  radius  vector  to  the  Earth  center 
If  the  total  estimator  velocity  estimate  exceeds  V  ,  there  is  a  high 
potential  for  generating  a  solution  which  can  skip  the  Earth 
atmosphere.  The  escape  velocity  is  the  J  2  times  the  minimum 
velocity  necessary  to  maintain  a  circular  orbit  at  the  given  altitude, 
or  radial  distance,  R.  A  decaying  Earth  satellite  has  generally 
followed  a  trajectory  which  has  spiralled  down  from  the  lowest 
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circular  orbit  which  could  be  maintained  in  the  presence  of  the 
atmospheric  drag.  (1)  Should  some  combination  of  the  scalar  fading 
memory  parameter,  y,  and  the  a  priori  state  covariance  matrix  yield 
a  solution  with  velocity  exceeding  the  circular  velocity  at  an  early 
epoch,  this  solution  may  be  discounted  on  physical  grounds  of  the 
true  system  dynamics.  While  it  does  represent  an  admissible  solution 
as  defined  by  the  magnitude  of  the  deweighted  state  covariance,  it 
will  generally  not  occur  within  the  true  reentry  dynamics.  One 
alternative  would  be  to  omit  these  initial  observations  from  the 
estimator  solution  and  begin  the  reentry  data  processing  later  during 
reentry . 

Fortunately,  the  Infinite  memory  formulation  shows  acceptable 
performance  in  these  regions.  Hence,  the  recommended  application  is  to 
restrict  the  use  of  fading  memory  until  later,  when  a  committed  to 
reentry  solution  is  the  only  admissible  estimator  solution.  This 
approach  still  allows  one  to  satisfy  the  basic  goal  of  carrying  the 
solution  forward  to  the  Earth  impact  locations  of  interest.  However, 
additional  research  is  recommended  to  develop  a  rigorous  examination 
of  the  stability  of  the  proposed  reentry  estimator  over  a  range  of 
reentry  conditions. 
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D.  Numerical  Simulation  Results 


D. 1 .  Deweighting  Technique  Variations 

The  series  of  single  sample  simulation  runs  summarized  in  Table 
XII  illustrate  some  relative  performance  considerations  of  a  number  of 
deweighting  methods  which  were  examined.  All  may  generally  be 
compared  to  the  basic  infinite  memory  estimator  results  to  discern  a 
gross  measure  of  the  estimator  performance  in  terms  of  error  and  the 
error/(ONE  SIGMA)  or  performance  ratio.  A  more  complete  indication  of 
the  estimator  performance  will  be  shown  later  with  the  Monte  Carlo 
results.  Complete  details  on  these  single  sample  results  are  contained 
in  Appendix  B.2. 

As  discussed  earlier  in  Section  A. ,  the  results  of  this  single 
sample  analysis  implementing  various  deweighting  methods  indicate  that 
the  adaptive  selection  of  a  scalar  deweighting  parameter,  y>  should 
offer  better  performance  than  the  infinite  memory  estimator,  a 
constant  scalar  deweighting,  or  an  adaptive  deweighting  of  only 
selected  states  within  S  , .  The  improvement  is  evident  in  both 
error  and  ratio  measures.  Only  in  the  constant  scalar  deweighting 
Case  2c,  does  performance  begin  to  compare  with  the  time  dependent 
adaptive  determination  of  y  at  each  epoch.  The  selection  of  such  a 
constant  deweighting  for  the  entire  trajectory  requires  some  a 
priori  knowledge  of  the  true  trajectory  dynamic  variations.  This  is 
partially  shown  by  comparison  of  Case  2b  and  2c  results.  A  y  =  1.56 
produced  better  results  than  the  y  =  1.10  of  Case  2b.  In  actual 
reentry  data  processing,  little  a  priori  information  for  proper 
selection  of  a  constant  y  will  exist. 
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Truth  Model:  8  =  f(Mach  no.)  +  1962  US  Standard  Atmosphere 
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Adaptive  Scalar  Deweighting 


techniques  on  the  digital  computer.  The  results  of  Table  XV  show  that  the  lea< 
coefficient,  a,  is  not  the  dominant  term  in  the  deweighting  scalar  expression. 


D.2.  Monte  Carlo  Results 
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With  the  implementation  of  an  adaptively  determined  scalar 
deweighting,  one  must  now  quantify  its  contribution  to  the  estimator 
performance.  A  sequence  of  Monte  Carlo  dynamics  mismatch  runs  and 
some  discussion  on  the  implementation  of  the  deweighting  scalar  will 
demonstrate  the  estimator  performance  in  terms  of  bias  and  RSS/(ONE 
SIGMA)  ratio.  A  Monte  Carlo  analysis  is  necessary  due  to  the 
complexity  of  the  nonlinear  dynamics  and  the  uncertainty  in  the  true 
reentry  dynamics  in  this  application.  These  problems  preclude  a 
simplified  analytic  development  of  bounds  on  the  estimator  performance. 
The  results  of  the  Monte  Carlo  analysis  lend  insight  into  the  viability 
of  the  proposed  technique  for  this  application. 

Three  cases  of  Monte  Carlo  results  were  completed  with  the 
fading  memory  estimator  formulation: 

1)  A  duplicate  of  the  infinite  memory  mismatched  dynamics  case 
now  with  an  adaptively  selected,  scalar  deweighting. 

These  runs  use  conditions  identical  to  the  infinite  memory, 
mismatched  dynamics  case  shown  in  Figures  14-17,  with  data 
from  a  single  observational  satellite  positioned  at  an 
initial  Right  Ascension  (RA)  of  +45°. 

2)  A  repeat  of  the  mismatched  dynamics  case  with  observational 
data  partially  overlapping  in  time  from  two  observation 
satellites  at  initial  R.A.  =  ±45°.  These  results  may  be 
compared  to  the  one  above  to  show  the  improvement  in 
estimator  performance  with  data  from  more  than  one 
observing  satellite. 


147 


3)  The  results  from  a  different  dynamic  mismatch  truth  model 
which  includes  a  step  change  in  3  late  in  the  trajectory. 

Data  comes  from  two  observational  satellites  initially 
positioned  at  ±45°  R.A.  This  case  was  selected  to  provide 
a  functional  representation  of  vehicle  fragmentation.  The 
dual  observations  occur  along  the  entire  data  span. 

In  all  cases,  the  same  baseline  as  used  in  the  infinite  memory 
Monte  Carlo  runs  apply.  This  includes  maintaining  the  same  initial 
seeds  on  the  state  variable  initial  condition  variations  and  on  the 
observation  noise  between  sets  of  the  30  Monte  Carlo  runs.  The  fading 
memory  results  will  contrast  the  estimator  performance  to  both  the 
"deweighted"  covariance  matrix  and  the  non-deweighted  covariance  matrix 
(Equations  111  and  112).  The  reference  information  for  the  three  sets 
of  Monte  Carlo  cases  is  shown  in  Table  XIII. 
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Table  XIII 


Fading  Memory  Monte  Carlo  Runs 


Truth 

Initial 

Observer 

Observer 
Data  Span 

Case  No. 

Model 

Trajectory 

R.A. 

(seconds) 

1)  IFAD 
(Single 
observer) 

ATM62 

8=f (Mach  no.) 

ALT  =  73.82KM 
R.A.  =  0° 

Decl  =  68° 
incl  I  10° 

+45° 

10  -  330* 

2)  2FAD 

ATM62 

ALT  =  73.82KM 

±45° 

1) 

40-300 

(Overlapping 

8=f (Mach  no.) 

R.A.  =  0° 

2) 

10-140 

dual 

Decl  =  68° 

observer  data) 

incl  2  10° 

3)  3FAD 

ATM62 

ALT  =  73.82 

±45° 

1) 

10-310 

(Dual 

observer 

data 

throughout) 

step  A8  at 

T  =  250 

R.A.  =  22° 

Decl  =  0° 
incl  ;  30.5 

2) 

10-310 

Observer  times  from  initial  epoch 
Observer  1)  at  R.A.  =  +  45° 
Observer  2)  at  R.A.  =  -  45° 
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The  fading  memory  selection  criteria  for  each  case  are  shown  in 


Table  XIV.  The  fading  memory  scalar  selection  criteria  differ  very 
slightly  from  case  to  case,  simply  to  explore  the  computer  mechaniza¬ 
tion.  As  the  later  numerical  results  will  show,  the  dominant 

contribution  is  from  the  max  (6x.//P  "  )  term  within  v»  and  not 

1  m,m-lii 

its  leading  numerical  coefficient,  a,  in  the  expression: 

y  =  (a  max  1 6x . I / /~P  "  )2  (120) 

1  1 1  m,m-l 

n 

The  form  of  the  expression  for  y  was  chosen  to  allow  the  estimator 
software  to  optionally  use  a  scalar  deweighting  or  the  more  general 
diagonal  deweighting  of  Equation  81. 

In  all  three  cases  the  same  relative  convergence  criterion  on 
the  residuals  is  applied.  From  Equation  34,  cj  =  .15.  This  relative 
convergence  criterion  was  empirically  determined  in  conjunction  with 
the  absolute  convergence  criterion  (a  from  Equation  35)  to  obtain 
small  changes  in  the  state  variable  estimates  on  the  final  iteration. 

The  fading  memory  was  not  allowed  to  operate  in  the  early,  near 
circular  orbit  trajectory  phases  to  prevent  divergence  in  the  initial 
estimator  solutions.  This  was  accomplished  by  maintaining  the 
infinite  memory  formulation  until  the  ratio  of  aerodynamic  to 
gravitational  acceleration  exceeded  a  certain  minimum  value,  as  shown. 

The  values  were  chosen  to  examine  the  sensitivities  of  the  estimator 
performance  in  the  upper  altitude  region.  With  the  use  of  the 
"escape  velocity"  test  discussed  in  the  stability  considerations  of 
Section  C,  a  systematic  means  to  select  this  value  will  be  available 
for  application  to  real  reentry  data. 
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Table  XIV 


Fading  Memory  Criteria 


CASE  IFAD 

1)  Do  not  apply  fading  memory  unless  the  ratio  of  aerodynamic/ 
gravitational  acceleration  >  0.10 

2)  The  absolute  convergence  criterion: 

mean  residual  <  0.5  x  one  sigma  observation  noise 

3)  Apply  fading  for  the  current  epoch  update  if : 

-  one  Sx.  >  1.05  /~P  I 

l  m,m-l . . 

Ii 


-  two  or  more  <5x.  >  /p  ' 

i  m,m-l . . 


ii 


4)  Fading  parameter,  y  =  { 1 . 1  [max  j  6x  .  |  / /p  .  ]}2 

i  m  •  rn*-  x  ,  , 

ii 

CASE  2FAD 

1)  Do  not  apply  fading  memory  unless  the  ratio  of  aerodynamic/ 
gravitational  acceleration  >  0.15 

2)  The  absolute  convergence  criterion: 

mean  residual  <0.5  one  sigma  observation  noise 

3)  Apply  fading  for  the  current  epoch  update  if: 


-  one  6x.  >  /  P " 

i  m.m-l^ 


4)  Fading  parameter,  y  =  {1.08[max|  <5x  ,  |  / T? 


m,m-l 


ii 
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CASE  3 FAD 


I 


1)  Do  not  apply  fading  memory  unless  the  ratio  of  aerodynamic/ 
gravitational  acceleration  >  0.5 

2)  Absolute  convergence  criterion 

mean  residual  <  0.6  one  sigma  observation  noise 

3)  Apply  fading  for  the  current  epoch  update  if : 

-  one  6x .  >  vP " 

l  ra,m-l . . 

n 

4)  Fading  parameter,  y  =  { 1 . 80 [max | 6x . | //  ]}2 

The  leading  coefficients  for  the  fading  memory  criteria  were  arbitrar¬ 
ily  selected  to  implement  the  estimator  on  the  digital  computer  without 
multiple  selections  of  a  scalar  fading  memory  at  any  given  epoch  point. 
Only  for  the  third  set  of  Monte  Carlo  results  was  the  value  of  the 
leading  coefficient,  a,  chosen  specifically  to  be  significantly 
greater  than  one  (a  =  1.80).  This  selection  was  made  to  illustrate 
the  potential  of  applying  sufficient  fading  such  that  the  ^  matrix 
could  yield  covariance  data  with  acceptable  RSS  performance. 


152 


The  results  of  the  Monte  Carlo  sets  are  shown  in  Figures  18-41. 
The  deweighted  covariance  one  sigma  values  were  used  for  plotting  the 
results  in  Figures  18-29.  For  comparison,  the  non-deweighted 
covariance  one  sigma  values  derived  from  Equation  112  are  shown  in 
Figures  30-41.  Because  of  the  dynamics  uncertainties,  the  deweighted 
variance  values  offer  a  more  conservative  estimate  of  the  uncertain¬ 
ties  in  the  state  solution  than  the  non-deweighted  expression.  In  all 
cases  the  non-deweighted  variance  data  shows  a  marked  improvement  in 
RSS  performance  when  compared  to  the  infinite  memory,  mismatched 
dynamics  results  of  Figures  14-17  in  Chapter  II.  However,  no  system¬ 
atic  technique  is  apparent  for  selecting  the  leading  coefficient,  a, 
such  that  a  zero-mean  bias  is  obtained  in  the  state  estimate.  With  a 
zero-mean  bias,  the  ^  expression  would  represent  a  minimum  variance 
covariance.  With  the  collection  of  data  from  many  reentries,  a 
systematic  technique  may  potentially  be  developed  for  selecting  a, 
such  that  a  zero-mean  estimator  solution  is  available.  This  poses 
much  the  same  problems  for  the  estimator  as  does  the  proper 
construction  of  a  pseudo-noise  compensation  to  the  dynamics.  The 

presentation  of  the  results  for  the  S  statistic  is  to  illustrate 

m,m 

the  potential  for  use  of  should  proper  tuning  be  successful. 

In  all  cases  the  results  for  the  estimator  are  displayed  in  the 
same  format  as  the  infinite  memory  Monte  Carlo  cases.  This  includes 
presentation  of  position  and  velocity  bias  and  variance  data  from  the 
estimator  solutions.  The  RSS  of  the  mean  square  errors  about  the  true 
solution  (Equation  75)  are  contained  in  the  figures  labeled 
"ESTIMATOR  PERFORMANCE"  within  the  RSS  term  of  the  ratio  RSS/ 

(ONE  SIGMA).  The  standard  deviation  about  the  mean  solution  (Equation 


76)  is  shown  only  in  Figures  20/21,  24/25,  28/29  to  illustrate  the 
consistency  with  the  estimator-computed  deweighted  covariance  data. 
Additional  Monte  Carlo  derived  standard  deviations  are  not  plotted 
for  ease  of  viewing  of  the  fading  memory  results.  In  both  cases  the 
standard  deviations  derived  from  these  Monte  Carlo  computed  covariance 
matrices  followed  a  trend  consistent  with  the  estimator-computed, 
deweighted  state  covariance  matrix.  However,  the  deweighted  state 
covariance  matrix  provides  a  conservative  estimate  of  the  RSS  error  in 
the  state  estimate.  The  use  of  the  standard  deviations  from  this 
deweighted  covariance  is  recommended  since  a  trajectory  dependent  bias 
still  exists  in  the  estimator  solutions. 

Figures  18-19  show  the  RSS/(0NE  SIGMA)  ratio.  As  before,  the  RSS 
is  derived  from  the  Monte  Carlo  data  (Equation  75).  The  ONE  SIGMA 
values  are  the  deweighted  covariance  position  or  velocity  uncertain¬ 
ties  obtained  from  the  application  of  the  fading  memory  estimator. 
Figures  20-21  show  the  estimator  mean  position  and  velocity  bias 
magnitudes  relative  to  the  average  value  of  the  estimator  deweighted 
covariance  "one  sigma"  or  standard  deviation  values. 

For  this  fading  memory  duplicate  of  the  infinite  memory 
mismatched  dynamics,  single  observation  satellite  case,  a  marked 
performance  improvement  is  evident,  compared  to  Figures  14-17.  Only 
very  near  to  the  last  observation  does  the  RSS/(0NE  SIGMA)  ratio  grow 
beyond  a  value  close  to  one.  VThen  this  does  occur  it  is  principally 
within  the  velocity  terms.  Similar  results  are  apparent  in  Figures 
20-21.  Only  in  the  velocity  term  does  the  bias  exceed  the  estimator 
one  sigma  levels,  and  then  on  just  two  observations.  By  the 
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final  observation,  the  bias  is  within  the  estimator  covariance  one 


sigma  level. 

In  Figures  20-21,  one  can  observe  the  onset  of  the  deweighting 
at  the  54  second  point,  and  the  general  ability  of  the  estimator  to 
minimize  the  divergent  behavior  of  the  infinite  memory  results  in 
Figures  16-17. 


ESTIMATOR  PERFORMANCE 


FIGURE  19 


FIGURE  10 
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Having  illustrated  an  improved  estimator  performance  for  the 
single  observation  satellite,  one  must  explore  the  higher  data  content 
and  improved  observability  offered  by  the  dual  observation  satellite 
data.  Figures  22-25  contain  these  results,  with  data  overlap  from  the 
two  observation  satellites  between  40  -  140  seconds  of  trajectory  time. 
This  data  overlap  region  was  selected  to  simulate  one  possible  orbital 
deployment  of  the  observation  satellites.  Doing  so,  it  includes 
regions  of  both  single  and  dual  observation  of  the  reentry  trajectory. 
As  the  results  in  Figures  22-25  show,  there  is  a  marked  improvement  in 
both  RSS/(0NE  SIGMA)  and  the  bias  magnitude  compared  to  the  deweighted 
covariance,  single  observer  results.  This  improved  performance  results 
with  the  dual  observations  of  the  trajectory  between  the  40  and  140 
second  points.  This  indicates  that  an  acceptable  estimator  solution 
should  then  be  available  for  propagation  to  Earth  impact.  The  bias  is 
well  below  the  standard  deviation  from  the  deweighted  covariance.  The 
deweighted  covariance  yields  a  conservative  measure  of  the  RSS  error 
for  the  position  and  the  velocity  terms. 
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FIGURE  22 
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The  final  set  of  Monte  Carlo  results  with  the  deweighted 
covariance  are  shown  in  Figures  26-29.  In  this  case,  a  different 
truth  model  was  used  to  generate  the  simulated  observations.  This 
included  a  step  change  by  50%  in  3  at  T  =  250  seconds  to  simulate 
reentry  object  fragmentation.  The  region  of  dual  observation  data 
included  the  entire  data  span  from  T  =  10  to  T  =  310  seconds.  The 
specification  on  fading  memory  selection  did  not  allow  fading  to  begin 
until  the  ratio  of  aerodynamic  to  gravitational  acceleration  was 
greater  than  0.5.  This  value  was  selected  to  allow  examination  of 
where  the  RSS/(0NE  SIGMA)  ratio  grew  beyond  a  value  of  one  with  the 
infinite  memory  formulation  for  this  different  truth  model  case. 

The  effect  of  this  restriction  is  evident  in  Figure  26-29. 

Prior  to  the  trajectory  time  of  approximately  174  seconds,  the  solution 
is  simply  an  infinite  memory  solution  with  dual  observational  data. 
After  174  seconds  of  trajectory  time,  the  adaptive  fading  memory 
selection  of  y  is  then  applied.  This  fact  is  shown  most  clearly  in 
Figures  26-27  where  a  slight  growth  in  the  RSS/(0NE  SIGMA)  ratio  is 
evident  prior  to  the  point  where  deweighting  commences.  After  this 
point,  the  RSS/(ONE  SIGMA)  ratio  shows  acceptable  performance  using 
the  deweighted  covariance  values.  The  RSS/(0NE  SIGMA)  ratio  is  less 
than  or  equal  to  one.  The  bias  in  position  and  velocity  is  well  below 
the  standard  deviation  from  the  deweighted  covariance.  The 
relatively  conservative  values  available  from  the  deweighted 
covariance  are  evident  in  Figure  27,  where  the  velocity  RSS/(0NE  SIGMA) 
ratios  are  consistently  much  lower  than  one  with  the  use  of  the 
fading  memory. 
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The  availability  of  dual  observations  through,  and  beyond,  the 
step  change  in  &  at  T  *  250  seconds,  enhances  the  ability  of  the 
estimator  to  accommodate  this  simulated  fragmentation  effect.  This 
verifies  the  speculation  raised  while  examining  the  infinite  memory 
single  sample  runs  of  Table  I.  That  is,  given  sufficient  data,  the 
estimator  can  adapt  its  local  model  to  large  discrete  changes  when 
employed  in  the  adaptive,  fading  memory  formulation. 
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Having  illustrated  the  validity  of  adaptive  selection  of  a 

scalar  deweighting  by  examining  the  bias  and  RSS  performance  relative 

to  the  deweighted  covariance  matrix,  one  must  also  consider  the 

magnitude  of  the  non-deweighted  covariance  matrix,  S  ,  and  its 

m,m 

ability  to  provide  acceptable  estimator  statistics.  Figures  30-41 
repeat  the  Monte  Carlo  results  for  the  non-deweighted  covariance 
matrix  (Equation  112).  One  can  see  a  marked  improvement  in  the  esti¬ 
mator  performance  with  the  non-deweighted  covariance: 

1)  Relative  to  the  infinite  memory  results,  and 

2)  Increasing  with  the  enhanced  observability  and  data  content 
of  the  dual  observer  cases  (Figures  34-41). 

In  the  first  two  cases,  the  RSS/ (ONE  SIGMA)  non-deweighted  ratio 

grows  slightly  above  one  near  the  end  of  the  data  span.  The  single 

observer,  mismatched  dynamics  case  has  several  data  points  where  the 

bias  exceeds  the  non-deweighted  covariance  one  sigma  (Figures  32  and 

33).  With  the  addition  of  dual  observational  data,  only  two  velocity 

bias  points  exceed  the  non-deweighted  covariance  one  sigma  levels 

(Figure  37).  Lastly  in  the  dual  data  case  with  a  step  change  in  g, 

one  can  see  the  potential  for  developing  a  deweighting  selection 

criterion  which  can  possibly  deliver  acceptable  estimator  statistics 

from  the  non-deweighted  covariance  matrix,  S  .  The  RSS/(0NE  SIGMA) 

m,m 

velocity  ratio  grows  beyond  a  value  of  one  (Figure  39) ,  while  the 
position  ratio  is  remarkedly  close  to  one  throughout  the  entire  data 
span  (Figure  38).  The  position  bias  is  consistently  inside  the  non- 
deweighted  covariance  one  sigma  magnitudes  (Figure  40).  The  velocity 
bias  only  exceed  the  non-deweighted  covariance  one  sigmas  at  the  very 
last  few  data  points. 
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One  must  recall  that  the  deweighting,  criterion  of  Table  XIV 

were  developed  solely  to  insure  acceptable  estimator  performance 

relative  to  the  deweighted  covariance  matrix.  Therefore,  all  the 

results  with  the  non-deweighted  covariance  matrix,  S  ,  shown  here 

m  ,m 

are  simply  a  by-product  of  those  scalar  deweighting  criterion.  The 

potential  does  exist,  however,  to  iteratively  select  the  scalar 

deweighting  with  a  more  stringent  selection  criterion,  such  as  the 

<5x,  <  v/~S  "  ,  or  6x.  <  bt,r_P  "  (b  ^  1.0).  Examination  of 

1  m,m-l . .  l  m,m-l . . 

n  n 

this  potential  was  not  completely  investigated  in  the  current  analysis, 
since  acceptable  estimator  statistics  were  available  from  the 
deweighted  covariance. 

One  must  also  recall  the  circumstances  under  which  the  non- 
deweighted  covariance  is  expected  to  yield  an  accurate  measure  of  the 
random  error  in  the  estimator  solution  for  the  state  trajectory.  In 
the  deterministic  dynamics  model  formulation,  the  model  must  be  exact 
to  yield  an  unbiased  estimate  of  the  trajectory  (13).  This  may  never 
be  the  case  when  applied  to  the  relatively  unknown  dynamics  of  real 
reentry  trajectories.  The  random  error  indicated  by  the  non-deweighted 
state  covariance  matrix  does  not  yield  a  good  definition  in  the  total 
estimator  error  due  to  the  bias  in  the  state  estimate. 

There  is  a  very  short  time  span  of  validity  for  the  deterministic 
dynamics  model  to  represent  the  truth  model  dynamics  locally  between 
epochs.  This  can  be  illustrated  by  considering  the  magnitude  of  the 
adaptively  selected  scalar  deweighting  parameter,  y.  One  may  define  a 
half-life  on  the  covariance  matrix  as  follows: 

^P  =  (l/y)AS  P  (121) 
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Solving  for  the  time  constant  yields: 


ln(S) 

ln(l/y) 


(122) 


This  it!  gives  an  indication  of  the  time  span  of  local  model  validity 
-2 

as  measured  by  the  half-life  of  the  a  priori  covariance  matrix  at  each 
epoch  where  a  new  y  is  adaptively  determined.  Table  XV  shows  typical 
estimator  runs  from  the  dual  observer  Monte  Carlo  runs.  The  magnitude 
of  y  and  the  approximate  half-life  from  Equation  122  are  listed. 

A  review  of  Table  XV  shows  the  following: 

1)  The  dominant  contribution  to  the  magnitude  of  the  adaptively 

determined  scalar,  y,  is  the  |<5x.|//p '  term:  and  not, 

‘  l  m  ,m-l  .  . 

n 

the  leading  numerical  coefficient,  a. 

2)  The  time  span  of  local  model  validity  is  very  short, 
indicating  one  should  not  anticipate  unbiased  estimator 
solutions. 

As  a  result  of  these  considerations,  one  would  not  recommend  the 

utilization  of  the  non-deweighted  covariance,  S  ,  as  an  indication  of 

m  ,m 

the  random  error  in  the  trajectory  solution.  In  light  of  the  uncertain 

dynamics  of  true  reentries,  a  systematic  means  to  modify  the  leading 

coefficient  within  the  expression  for  y  to  drive  the  RSS  performance 

within  the  standard  deviation  of  S  is  also  not  apparent.  Since 

Tn  ,ra 

acceptable  RSS  performance  was  evident  with  the  deweighted  covariance 
matrix,  it  is  considered  prudent  to  use  the  estimator  deweighted 
covariance  matrix  for  a  measure  of  the  uncertainty  in  the  estimator 


solution. 


Table  XV 


Fading  Memory,  Covariance  Half-Life 


CASE  2FAD  B  =  f (M)  +  ATM  62 


Y  =(1 . 08  max  j  <5x^  | 

//  p  ,  )2 
m,m-l . . 
n 

Approximate 

Covariance 

Time 

Y 

Half-Life 

84. 

2.907 

.650  (secs' 

94. 

5.503 

.406 

104. 

2.634 

.710 

124. 

1.935 

1.050 

134. 

1.545 

1.590 

154. 

1.633 

1.413 

214. 

3.  920 

.507 

254. 

5.253 

.418 

274. 

2.409 

.788 

284. 

3.717 

.528 

294. 

7.161 

.352 

CASE  3FAD  AB  at  T  =  250  +  ATM  62 


Y  =  (1.80  max  I  <5x  .  | //~P  "  )2 

x  m,m-i . . 


ii 

164. 

33  .  907 

.197 

174. 

21.455 

.226 

184. 

15.179 

.254 

204. 

14. 296 

.261 

214. 

16.704 

.246 

234. 

20.784 

.229 

254. 

26.420 

.212 

284. 

48.066 

.180 

294. 

22.734 

.222 

304. 

44.676 

.182 

Note: 

omitted  epoch  points  at  10 

sec  data  rate  did  not 

generate  additional  fading 

,  via  the  dynamically 

selected 
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The  deweighted  covariance,  P  ,  has  characteristics  which  are 

m  ,m 

analogous  to  the  covariance  from  a  properly  tuned  estimator,  with 
pseudo-noise  strength  chosen  such  that  the  estimator-computed  and  true 
mean  square  errors  agree  well. 

These  considerations  are  particularly  apparent  when  one 
considers  the  magnitude  of  the  ultimate  covariance  matrix  when  propa¬ 
gated  to  impact.  The  results  of  the  next  section  show  a  marked  reduc¬ 
tion  in  even  the  deweighted  covariance  matrix  one  sigma  position 
uncertainties  relative  to  previous  operational  experience.  The 
deweighted  covariance  provides  an  acceptable  and  conservative  measure 
of  the  uncertainties  in  the  ultimate  satellite  impact  point. 


I 


D. 3.  Propagation  to  Impact 

Table  XVI  shows  various  methods  of  propagating  the  two  mismatched 

dynamics  cases  using  dual  observation  data  to  impact  with  one  sigma 

uncertainties  for  each  state  variable.  Similar  to  the  infinite 

memory  impact  covariance  results,  the  linear  propagation  methods  do 

not  reflect  a  completely  accurate  portrayal  of  the  second  order 

statistics  at  the  ultimate  impact  point.  In  both  cases,  the  non- 

deweighted  variance  data  have  one  sigma  values  which  are  smaller  than 

the  average  bias  in  selected  state  variables.  This  is  particularly 

true  for  the  velocity  terms.  It  is  also  true  for  the  Bp  term, 

o 

principally  since  the  CASE  2FAD  truth  model  had  a  variable  gp^.  Also, 
the  discrete  step  change  of  6pq  in  CASE  3FAD  still  reflects  the  mis¬ 
match  between  the  estimator  and  the  truth  model. 

With  the  selected  deweighting  criteria,  the  accuracy  of  the  non- 
deweighted  covariance  values  is  anticipated  to  be  a  less  conservative 
measure  of  the  random  state  uncertainty  than  that  available  from  the 
deweighted  state  covariance.  The  linear  propagation  of  the  deweighted 
state  covariance  is  often  overly  conservative,  with  very  large  one 
sigma  values  relative  to  the  bias  on  selected  state  variables.  This  is 
particularly  true  for  the  final  case  which  included  a  step  change  in  6, 
where  the  deweighted  one  sigmas  are  much  larger  than  the  bias  solution 
for  all  state  variables.  For  these  reasons,  a  Monte  Carlo  derived 
impact  covariance  for  the  final  propagation  phase  is  recommended.  The 
ratio  of  the  RSS  of  the  mean  square  errors  to  one  sigma  about  mean 
solutions  provide  the  best  indication  that  a  Monte  Carlo  derived 
covariance  provides  the  most  accurate  reflection  of  the  RSS  uncertainty 


in  the  state  solution  at  Earth  impact. 


Fading  Memory  Estimator  -  Impact  Statistics 
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In  both  these  cases,  the  impact  uncertainties  are  orders  of 
magnitude  less  than  previous  operational  experience  of  propagating  the 
last  orbital  estimate  to  impact.  This  offers  the  most  convincing 
evidence  of  the  chosen  deweighting  estimation  technique  affording 
significant  improvements  to  previous  impact  location  and  uncertainty 
determinations. 

Figures  42  and  43  show  the  Earth  tangent  plan  projections  of  the 
estimator  solutions  for  two  cases  of  dynamics  mismatch  relative  to  the 
estimator  dynamics  model.  The  mean  estimator  impact  location  is  shown 
as  point  0,  with  the  truth  model  impact  at  point  X.  The  difference  in 
distance  between  these  points  show  the  bias  in  the  ultimate  impact 
solution.  Both  one  (lo)  and  two  (2o)  standard  deviation  error  ellipses 
from  the  "deweighted"  state  covariance  matrix  are  plotted.  These  were 
developed  by  a  Monte  Carlo  propagation  of  trajectories  from  the  final 
reentry  data  point.  As  can  be  seen,  consistent  impact  locations  and 
statistical  descriptions  of  the  uncertainties  are  provided.  The 
magnitude  of  bias  is  less  than  the  one  standard  deviation  magnitude. 

The  magnitude  of  the  position  uncertainties  is  significantly  improved 
over  the  many  hundreds  or  thousands  of  kilometer  uncertainty  of  Earth 
impact  available  from  existing  operational  methods  of  such  agencies  as 
the  USAF  SPACETRACK  System. 

Figure  44  provides  an  illustration  of  how  large  the  impact 
uncertainties  would  be  with  a  simple  propagation  of  the  initial  epoch 
(73.82  km  altitude)  covariance  to  impact.  The  initial  epoch  contained 
a  standard  deviation  of  position  of  approximately  3.5  km,  mostly  in¬ 
track  error  along  the  trajectory.  This  results  in  over  a  400  km  one 
sigma  uncertainty  along  the  path  of  the  reentry  ground  trace,  at 
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FIGURE  44 


Earth  impact.  Figure  44  may  be  directly  compared  to  the  one  sigma 
error  ellipse  of  Figure  42  to  see  the  significance  of  the  proposed 
technique  in  identifying  the  satellite  debris  search  area. 
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Chapter  IV  -  Cone  I  unions 

ihe  principal  goal  of  this  research  was  not  to  enhance  the 
existing  estimation  theory  or  reentry  dynamics  theory  by  a  significant 
amount .  The  principal  goal  was  to  bridge  the  considerable  gap  between 
theory  and  engineering  practice.  A  valuable  tool  for  astrodynami.es 
research  has  been  developed  as  a  product  of  this  dissertation.  One 
could  have  foresaken  this  effort  by  saying  that  the  original  research 

constraints  were  too  severe.  The  use  of  an  orbital  observer  as". . 'p.d 

angles  only  (IR)  data  at  a  10  second  data  rate.  The  angular 
observations  provide  a  limited  observability  of  the  system.  The  10 
second  data  rate  does  not  capture  many  of  the  rapidly  varying  dynamics 
changes  of  reentry.  The  research  approach  concentrated  on  extending 
existing  orbit  estimation  methods,  which  often  use  a  deterministic 
dynamics  model,  into  the  reentry  regions.  The  deterministic  estimator 
formulation  is  difficult  to  apply  in  an  uncertain  dynamics  region. 
Perhaps  this  is  why  this  particular  problem  has  perplexed 
astronautical  engineers  for  the  better  part  of  20  years.  With  the 
collection  of  many  additional  true  reentry  trajectory  estimates,  the 
issues  presented  in  this  dissertation  can  be  married  to  more 
sophisticated  estimation  theory  to  provide  further  advances  in  this 
complex  area  of  research. 

A  basic  objective  or  this  research  was  to  extend  the  existing 
orbit  determination  methods  into  the  high  dynamic  model  uncertainty 
regions  of  reentry.  A  differential  correction  estimation  technique 
was  utilized  to  process  angular  infrared  (IR)  observations  of  a 
reentry  trajectory  to  reconstruct  the  Larth  impact  locations  and 
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uncertainties  of  arbitrary  decayed  satellite  trajectories.  The  azimuth 
and  elevation  1R  observations  were  obtained  from  an  orbiting  sensor 
viewing  the  reentry  from  a  synchronous  orbit.  The  linearized, 
differential  corrector  used  a  deterministic  dynamics  model  in  a  Taylor's 
series  expansion  about  a  reference  trajectory.  The  technique  provided 
state  variable  updates  and  covariance  data  for  the  reentry  trajectory 
positions,  velocities,  a  vehicle  ballistic  parameter,  and  the  density 
scale  height  of  an  exponential  atmospheric  density  model. 

Previous  operational  experience  (5)  has  shown  that  extreme 
difficulty  exists  with  application  of  a  pseudo-noise  compensation  to 
the  estimator  model  dynamics  for  orbital  satellite  applications.  Up  to 
100  Earth  revolutions  of  tracking  data  of  the  orbiting  satellite  have 
been  required  to  tune  the  dynamic  pseudo-noise  coefficients  and 
covariance  matrix  elements  properly  for  a  specific  satellite  of 
interest.  With  the  short  time  spans  and  limited  tracking  data  of  the 
reentry  trajectories,  this  task  is  even  more  difficult  to  implement. 
Using  a  fading  memory,  recursive  estimator  formulation,  the  adaptive 
estimation  of  an  ad  hoc  scalar  fading  memory  parameter  provided 
estimates  significantly  improved  over  existing  techniques  for  reentry 
trajectory  and  Earth  impact  estimation. 

The  structure  of  this  technique  is  similar  to  using  a  scalar 

tuning  parameter  to  adaptively  tune  a  noise  covariance  such  as 

Q(t)  =  Q  (t)  +  Yt\Q(t).  In  the  current  formulation,  the  estimator- 
o 

computed,  or  a  priori,  covariance  is  utilized  as  the  starting 
definition  of  uncertainty,  and  modified  by  the  adaptively  determined 
scalar  parameter,  y. 
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Specific  conclusions  of  the  research  included  the  following: 

1)  Monte  Carlo  simulation  analyses  show  that  the  dynamics 
uncertainties  of  the  general  satellite  decay  trajectories  to 
significantly  affect  the  estimator  performance.  Without 
model  compensation,  significant  bias  exists  in  the  state 
variable  solutions  relative  to  the  standard  deviation  of  the 
estimator-computed  state  covariance  matrix. 

2)  With  a  deterministic  dynamics  model,  angular  observation 
accuracies  with  standard  deviations  less  than  10-5  radians 
induce  significant  error  in  the  state  estimate  relative  to 
the  standard  deviations  of  the  state  covariance  matrix. 

3)  Due  to  the  dynamics  uncertainties  anticipated  in  real  reentry 
applications,  a  recursive  formulation  of  the  estimator  is 
recommended  which  uses  a  short  time  span  between  the 
trajectory  update  point  (epoch)  and  the  observation (s)  being 
processed.  Monte  Carlo  analyses  of  the  estimator  bias  and 
RSS  performance  demonstrated  that  successful  linearization 
(as  defined  by  Equation  44)  relative  to  the  reference 
trajectory  is  obtained  by  limiting  the  time  span  between 
epoch  update  point  and  the  observations  used  for  this 
trajectory  update. 

4)  For  acceptable  estimator  performance  in  terms  of  bias  and 
RSS/ (ONE  SIGMA)  ratio,  observations  from  more  than  a  single 
orbital  source  are  required  to  provide  higher  data  content 
over  similar  time  spans  and  improved  observability  of  the 
reentry  trace.  The  additional  use  of  range,  or  range  rate, 
observations  should  improve  the  observability  of  the  system. 
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5)  An  eight  dimensional  formulation  (3  positions,  3  velocities, 
ballistic  parameter,  density  scale  height)  provided  superior 
estimator  performance  to  a  seven  dimensional  formulation. 

The  standard  orbit  determination  approach  which  uses  a  seven 
dimensional  state  vector  (3  positions,  3  velocities, 
ballistic  parameter),  with  density  from  a  standard  density 
model,  proved  inadequate  in  the  current  estimator  formulation. 
The  improvement  in  the  performance  for  the  eight  dimensional 
system  derives  from  simpler  mathematics  in  the  estimator 
dynamics  model  and  continuous  valued  partial  derivatives  of 
the  dynamics  over  the  trajectory  space  for  the  Taylor's 
series  linearization. 

6)  With  a  deterministic  dynamics  model,  the  time  span  of  model 
validity  relative  to  true  reentry  dynamics  is  very  short.  A 
very  limited  number  of  observations  are  available  over  this 
short  time  span.  A  proper  statistical  testing  of  these 
observation  residuals  is  not  available  for  adaptive 
compensation  of  the  estimator  dynamics.  A  magnitude  check 
on  the  state  variable  changes,  <5x^,  at  each  update  epoch 
relative  to  the  standard  deviation  of  the  state  covariance 
matrix  provides  a  satisfactory  measure  to  determine  which 
trajectory  points  required  additional  fading  on  the 
estimator  memory.  This  is  demonstrated  by  examining  the 
bias  and  RSS  performance  of  the  estimator  in  a  Monte  Carlo 
analysis.  This  adaptive  selection  of  a  scalar  fading  memory 
parameter  is  easily  incorporated  into  the  estimator  structure. 
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7)  Lastly,  a  Monte  Carlo  propagation  to  impact  preserves  the 


integrity  solution  statistics  over  the  final  non-observable 
portion  of  the  trajectory.  The  Monte  Carlo  derived 
covariance  about  the  mean  and  the  true  impact  solutions 
compare  very  closely.  Both  these  measures  yield  values 
significantly  different  than  those  obtained  from  a  linear 
propagation  of  the  final  epoch  covariance  to  impact. 

Application  to  general  decay  trajectories  of  the  fading  memory 
estimator  provides  impact  locations  and  covariance  data  consistent  with 
the  RSS  uncertainties  of  those  locations.  The  magnitudes  of  the  impact 
location  uncertainties  are  between  one  to  two  orders  of  magnitude 
smaller  than  those  available  from  current  operation  techniques.  The 
method  is  applicable  to  world-wide  impact  locations  with  the 
availability  of  the  orbital  reentry  observations. 

Users  of  the  proposed  estimation  algorithm  are  cautioned  to  refer 
to  the  stability  considerations  of  Chapter  III,  Section  C.  A  rigorous 
stability  analysis  has  not  been  completed  for  the  estimator  algorithm. 


Chapter  V  -  Recommendations 


1)  Investigations  which  vary  the  observation  data  rate  and  the 
time  varying  character  of  true  reentry  dynamics  are 
necessary  to  examine  the  estimator  performance  extensions. 
This  should  be  conducted  with  an  investigation  of  using  more 
accurate  observations  (much  less  than  10~5  radians)  in 
conjunction  with  higher  data  rates  and  higher  frequency 
variations  in  reentry  dynamics.  Examination  of  using 
alternative  measurements,  such  as  range  or  range  rate,  is 
also  recommended. 

2)  Further  analysis  of  a  means  to  apply  the  fading  memory 
estimator  to  very  high  altitude  and  shallow  reentry  angle 
reentries  is  recommended.  Should  violent  dynamics  changes, 
such  as  vehicle  fragmentation,  occur  in  these  very  early 
data  regions,  the  suggested  fading  memory  scalar 
determination  may  result  in  divergent  estimator  performance 
due  to  atmospheric  "skip"  trajectories  being  admissible 
solutions  within  a  large  magnitude  deweighted  state 
covariance  matrix.  A  rigorous  stability  analysis  approach 
needs  to  be  defined  for  both  the  state  update  and 
propagation  equations. 

3)  Further  investigation  of  implementing  the  fading  memory 
would  be  of  interest.  It  is  desireable  to  minimize  the 
bias  in  the  estimator  solution  such  that  the  non-deweighted 
state  covariance  matrix  provides  a  minimum  variance 
covariance  for  the  state  errors  (with  zero  mean  <5x(t) 
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estimates).  However,  since  the  size  of  the  impact 
uncertainties  from  the  deweighted  state  covariance  matrix 
provide  such  a  marked  improvement  over  existing  methods, 
they  are  acceptable  for  most  reentry  debris  search  and 
recovery  requirements. 

4)  The  application  of  this  suggested  estimation  technique  to  a 
wide  class  of  reentry  trajectories  should  provide  a  large 
empirical  data  base  for  improvement  of  the  estimator 
dynamics  model.  Subsequent  investigations  should  then  be 
pursued  such  as  dynamic  model  pseudo-noise  compensation, 
statistical  linearization,  or  higher  order  filters.  These 
may  then  be  compared  to  the  recommended  fading  memory  tech¬ 
nique  for  comparison  of  relative  performance. 
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Appendix  A  Partial  Derivatives 


A. 1.1.  Partial  Derivatives  of  Drag  Accelerations 


Recall  d^,  ,  d^  =  drag  accelerations  in  x,  y,  z  components, 

therefore  take  the  partial  derivatives  by  chain  rule  - 

3/3x.  H*/Q: 

1 


i;(H*/Q)  -  GM/qgo 


R  (1-f)  (2xR2 (l-2f+f 2 )  +  xz2 (2f-f 2 ) 
o _ _ 

R3 ( (l-2f+f 2) (x2  +  y2)  +  z2)3/z 


a  *  R  d-f)  (2yR2 (l-2f+f ‘ )  +  yz2(2f-f2)  y 

ar(H  /Q)  =  GM/Qg  — 5 - 7 - 

3y  °  R3((l-2f+f2)(x2  +  y2)  +  z2)3/2  R3 


(A4) 


3  * 

li  (H  /Q)  " 

R  (1-f)  (2zR2 (l-2f+f 2)  +  zR2 (2f-f 2 )  +  z3(2f-f2)  z 

GM/Qg  - - - y - 

R3 ((l-2f+f2) (x2  +  y2)  +  z2)3'2  R3 

VR  =  [(x  +  ujy)2  +  (y  -  uxf  +  z2]'5  (A5) 

3/3x.  VR: 

3VR/3x  *  -u(y-u'x)/VR 
3Vr/3x  =  (x  +  uiy)/VR 

3VR/3y  =  oj(x  +  tuy)/VR  (A6) 

3VR/3y  =  (y-wx)/VR 

3V  /3z  =  0 
K 

3Vr/3z  *  z/VR 
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Drag  Partials: 


3fd 

=  -^6P  (x+my )  [V  3/3x  e"H  /Q  +  e' H  /Q  3/3x  VI 

dX  O  K  K 


(A7 ) 


3f 


d  *  , 

x  i  q  -H  /Q 
—  =  -4 Bp  e  V 


3x 


•k 

-^Bp  e  H  ^  (x+-y)  3/3x  V 

O  K 


(A8) 


3f 


ux  lB  -H*/Q  v 
3y~  =  '4Bpo  6  VR 

*  * 

-kBp  (x+uy)  [V  3/3y  e' H  /Q  +  e“H  /Q  3/3y  V  ] 

OK  K 


(A9) 


3f 


3y 


k 

—  =  -hBp  e"H  ^  (x+uy)  3/3y  V 

O  K 


(A10) 


3fa  * 

•  -  H  /0 

— jjj-  =  -^BPp  (x+wy)  VR  3/3z  e 


(All) 


3f 


9z 


* 

—  =  -*s6p  e  H  ^  (x+uy)  3/3z  V 

•  O  K 


(A12) 


3f 


d  * 

x  i  -h  /q  ,,  ,•  v 

—  =  -H  e  VR  (x+wy) 


3(BP0) 


(A13) 


3f 


X  _  _-H*/Q 


3(Q)  -  -^o  VR  <^>  H  ^ 


(A14) 


3f 


r-2  -  SfBp  e" H  /Q  V  u> 
3x  o  R 


(A15) 


3/3x  e'H  /Q  +  e‘H  /Q  3/3x  VD] 


-h&PQ  (y-wx)  [ VR  ^  —  -RJ 
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(A16 ) 


1  d  * 

— ^  =  -JsBp  e"H  /Q  (y-ojx)  3/  3x  V. 
«  *  0  J 
3x 


-1i6po  (y-ux)[VR  3/  3y  e 


-  H*/Q 
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+  e‘H  /Q  3/3y  VR] 


3f 
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-X  ,  J,«„  /« 
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3fd  *  * 


=  -^Sp  (y-uix)[V  3/3z  e  H  ^  +  e  ^  ^  3/3z  V  ] 
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3f 


3z 


d  /  -V 

-X  =  e  ^  (y-ux)  3/3z  V, 


R 


3f 


3(6P 
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3f 


3(Q) 
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i  =  _4Bp  e-H  /Q  (y_wx)  V  H/Q2 
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•  • 
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(A23) 
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(A26) 


(A27) 


(A28) 


(A29) 


(A30) 
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A. 2. 2.  Geopotential  Partial  Derivatives:  (Ref:  DMAAC  Network  Analysis 
Program,  Vol  I,  Part  I  -  Mathematical  Analysis,  Sep  1970) 

Define  partial  derivative  matrix  of  geopotential  function  W.R.T. 
state  variable  position  coordinates: 


) 

<P 

XX 

xy 

XZ 

yx 

<p 

yy 

$ 

yz 

! 

<P 

ZX 

zy 

zz 

where  <J>  =  34>/Sx3y  etc.  (A31) 


Define  ujA-jj  where  j  =  rotation  matrix  back  to  E.C.I.  Greenwich 
Meridian  reference  (note:  actually  required  only  if  tesseral  harmonics 
are  used) 


<p  =  1/4R  2  I  {C  [U®t;  +  cm  (C° J  Um"^  -  2Um, ,)]  + 
xx  o  n,m  n,m  n+2  n  n+1  n+2  n+2 


S  +  cm  (C®-1  vm-2  _  2V”\,)]} 

n,m  n+2  n  n+1  n+2  n+2 


t>  =  l/R  2  I  Cm  [C  U“  +  S  V™,] 
zz  o  n  n,m  n+2  n,m  n+2 

n  ,m 


(<fr  +  <t>  ) 

XX  zz 


P  =  <P  =  1/2R  2  l  -  (n-m+1)  [C  (C™  Um*l  -  U®**)  + 
xz  zx  o  n,m  n+1  n+2  n+2 

n  ,m 


S  (Cm+1  +  v®rb] 

n,m  n+1  n+2  n+2 


?)Um 
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1  -  (-2 (n-m+1)  uJJ+1) 
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3V 
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2R 


where:  C™  =  (n-ro+1)  (n-tn+2) 
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Geometry  Partial s 


i. 

{ 

t 


l 


A. 2. 


hj  =  sin-1  y 1  +  r. 

(y’-'+z'2)1' 

h;  =  sin-1  x  ' _ 

(x ' 2+y  '"+z  ,2)^ 

(Prime  ?  Observer  Coordinate  Frame) 


By  simply  chain  rule  H  = 

hn  =  3hi/3x  =  3hi/3x'  3x'/3x  + 

hi3  =  3h  i  /  3y  =  3hi/3x'  3x'/3y  + 

h  -  5  =  3hi/3z  =  3hi/3x'  3x'/3z  + 

h2;  =  3h;/3x  =  3h2/3x'  3x'/3x  + 

h23  =  3h2/3y  =  3h2/3x’  3x'/3y  + 

h2s  =  3h2/3z  =  3h2/3xf  3x'/3z  + 

All  other  his  elements  = 


3h/3x_^ 

3hi/3y'  3y'/3x  +  3hi/3z'  3z'/3x 

3hi/3y'  3y'/3y  +  3h:/3z'  3z'/3y 
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Numerical  Data 
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B.1.3.  High  Density  Exponential  Multipliers 
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B. 2.  Single  Sample  Analysis  Results 
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SINGLE  SAMPLE  ANALYSIS 
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INFINITE  MEMORY,  DYNAMIC  MISMATCH  (TABLE  I) 
TRUTH  MODEL:  •  ASSUMPTIONS 
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DEWEIGHTING  STUDIES  (TABLE  XII ) 

•  TROTH  MODEL!  •  ASSUMPTIONS:  ADAPTIVE  SCALAR  DEWEIGHTING 


SAP  III  Earth  Model  Parameters 

RE  =  6378.14  KM  Mean  Equatorial  Radius 

g  =  9.798  X  10"3  KM/sec2 
o 

f  =  1/229.256  Flattening  Factor 

GM  =  398601.3  Gravitational  Constant 

1 .  Zonal  and  Tesseral  Harmonic  Coefficients 
Normalization  Factor  =  -  /  2n+l 
Where  n  =  coefficient  subscript,  i.e.,  C^  or 
Example  Term: 

Unnormalized  ^2  =  ^2  0  = 

Normalization  Factor  =  - 
Therefore,  J 2  =  2.1653  X 


-  4.8417  X  10" 4 
/4+1  =  -  fT~ 

10- 4 
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